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Abstract 


This  dissertation  is  a  documented  case  history  in 
applying  an  analytical  method  to  a  field  problem.  First,  the 
correct  usage  of  the  existing  computer  program  ADINA 
(Automatic  Dynamic  Incremental  Nonlinear  Analysis)  is 
presented.  The  application  of  ADINA  to  landslide  problems  is 
explored . 

Strain-softening  stress-strain  behaviour  and  the 
resulting  progressive  type  failure  are  common  in  many 
landslide  problems,  but  the  analyses  available  are  empirical 
in  nature.  For  this  analysis,  a  load-transfer  technique  is 
developed  for  understanding  the  progressive  failure  process. 
The  material  model  used  results  from  a  strain-weakening 
approach. 

Examples  are  used  to  verify  the  load-transfer  program. 
This  approach  is  applied  to  a  simple  model  to  study  the 
stress  variation  along  the  slip-surface  (shear  zone). 

Finally,  the  approach  is  applied  to  the  Edgerton  Slide. 
The  numerical  results  from  the  case  study  are  compared  with 
the  available  laboratory  results  and  the  field  measurements. 
The  Young’s  Modulus,  used  in  the  calculating  displacements 
which  are  comparable  to  the  field  measurements,  agrees  with 
the  value  obtained  from  laboratory  procedure. 

The  applications  of  this  technique  to  other  areas  are 
discussed.  The  load-transfer  program  can  apply  to  other 
material  models  such  as  strain-stiffening,  creep  and 
no-tension  behaviour. 
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1.  Introduction  :  Scope  Of  This  Study 


1 . 1  Purpose  Of  The  Research 

The  basic  aim  is  to  apply  the  analytical  technique  to 
gain  an  insight  into  landslide  problems.  The  proposed 
technique  is  applied  to  an  analysis  of  the  Edgerton  Slide. 
This  case  history  has  been  studied  by  Robin  Tweedie  (1976) 
and  Ronald  Mokracki  (1982).  The  present  research  program 
involves  the  use  of  the  general  purpose  linear  and  nonlinear 
finite  element  analysis  computer  program  ADINA  (Automatic 
Dynamic  Incremental  Nonlinear  Analysis)  as  a  basic  tool.  The 
result  will  be  a  documented  case  history  applying  the 
analytical  method  to  a  field  problem. 

The  analytical  method  is  used  to  achieve  the  following 
goals : 

1.  Development  of  a  suitable  material  model  for  the  study 
of  the  shear  zone  behavior; 

2.  Development  of  an  analytical  procedure  for  the  study  of 
the  Edgerton  Slide; 

3.  Assessment  of  the  analytical  results  (primarily 
comparing  field  data  with  the  computer  output); 

4.  Assessment  of  the  general  application  of  the  proposed 
technique  (i.e.  load-transfer  technique). 
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1.2  Organization  Of  The  Thesis 

The  remaining  part  of  this  Chapter  will  deal  with  a 
literature  review  of  the  nature  of  progressive  failure  and 
some  analytical  work  dealing  with  progressive  failure. 

Chapter  2  presents  a  description  of  the  site  and 
updated  field  data  (to  August , 1 982 ) . 

The  formulation  of  the  analytical  methods  is  given  in 
Chapter  3,  primarily  to  explore  and  understand  the  available 
functions  of  ADINA.  The  procedures  developed  will  couple 
ADINA  with  the  load-transfer  program. 

Chapter  4  presents  the  Finite  Element  analysis  of  the 
Edgerton  Slide. 

Chapter  5  offers  the  general  application  of  the 
load-transfer  technique  and  the  conclusions  of  this  thesis. 

1.3  Review  Of  Nature  Of  Progressive  Failure 

Terzaghi  (1936)  suggested  that  the  removal  of  lateral 
support  in  stiff  fissured  clays  could  cause  opening  of  the 
fissures.  Moisture  ingress  leads  to  a  reduction  in  average 
strength  and  allows  more  deformation.  This  was  supported  by 
Cassel  (1948)  and  Skempton  (1948). 

Terzaghi  and  Peck  (1948)  and  Taylor  (1948)  pointed  out 
that  non-uniform  straining  of  strain-softening  material 
cannot  obtain  full  peak  strength.  The  soil  along  part  of  the 
sliding  surface  may  be  exerting  its  peak  strength  while  that 
along  the  remainder  may  be  exerting  a  smaller  value.  This 
hypothesis  forms  the  basis  of  the  definition  of  progressive 
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failure . 

Skempton  (1964)  in  the  Fourth  Rankine  Lecture 
introduced  the  phenomenon  of  residual  strength  which  leads 
to  the  question  in  slope  analyses  :  what  strength  parameters 
(i.e.  peak  strength  or  residual  strength  )  should  be  used  in 
the  design  of  slopes? 

Bjerrum  (1967)  in  the  Third  Terzaghi  Lecture  postulated 
a  mechanism  for  progressive  failure  as  a  result  of  a  large 
content  of  "recoverable  strain  energy"  in  overconsolidated, 
plastic  clays.  The  conditions  for  this  mechanism  to  occur 
were  : 

1.  The  material  shows  a  large  and  rapid  strength  decrease 
after  maximum  strength  is  exceeded. 

2.  Local  shear  stresses  tend  to  exceed  the  maximum 
strength . 

3.  Large  movement  due  to  the  release  of  locked  in  strain 
energy . 

Bishop  (1967)  suggested  a  mechanism  based  on  local 
overstressing  in  terms  of  the  shear  stress  (undrained 
condition)  or  the  ratio  of  shear  stress  to  the  effective 
normal  stress  (drained  condition).  After  the  formation  of  a 
zone  of  plastic  equilibrium  at  one  point  in  the  slope,  the 
zone  of  failure  would  propagate  along  the  potential  slip 
surface . 

Skempton  and  Petley  (1967)  suggested  that  large  strains 
would  require  a  pre-shearing  of  the  clay  and  the  forming  of 
principal  shear  planes.  Such  movement  would  be  obtained  by 
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processes  such  as  landsl iding ,  tectonic  movements  and 
glacial  ice  movements.  Meantime,  Peck  (1967)  pointed  out 
that  a  major  factor  affecting  our  ability  to  predict  whether 
or  not  a  slide  will  occur  is  whether  or  not  we  are  in  an  old 
slide  area. 

Yudhbir  (1969)  concluded  that  the  release  of  horizontal 
stress  (Ko  effect)  in  these  soils  is  a  dominant  factor. 
Bishop  and  Lovenbury  (1969)  showed  that  long  term  loading 
does  not  lead  to  substantial  strength  reductions,  which 
suggests  that  there  is  no  path  to  residual  which  by-passes 
the  peak. 

James  (1971)  pointed  out  that  large  displacements, 
often  in  the  order  of  feet,  are  necessary  to  develop 
residual  conditions  on  a  continuous  slip  surface. 

Morgenstern  (1977)  pointed  out  that  there  appear  to  be 
no  well-documented  case  histories  of  first  time  slides  in 
heavily  overconsolidated  soils  to  indicate  that  progressive 
failure  plays  a  dominant  role  in  governing  stability. 

The  preceding  discussion  indicates- that  stiff-fissured 
clays  and  clayshales  present  difficult  slope  stability 
problems.  These  difficulties  can  be  summarized  as  the 
stress-strain  relationship  of  the  soil,  the  effects  of 
fissures  and  openings  and  the  large  horizontal  stresses. 
However,  case  histories  show  that  many  slope  failures  in 
stiff  fissured  clays  and  clayshales  cannot  be  explained  in 
terms  of  peak  strength  values  and  equilibrium  methods  of 
analysis.  Therefore,  finite  element  method  may  help  us  to 
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gain  an  insight  on  the  failure  mechanism  of  stiff  fissured 
clay  slopes. 

1.4  Review  Of  Analytical  Work  In  Modelling  The  Problem  Of 
Progressive  Failure 

Prior  to  1970,  most  of  the  analyses  were  based  on  limit 
equilibrium  methods  such  as  the  Simplified  Bishop,  the 
Simplified  Janbu,  and  the  Morgenstern  and  Price  methods. 

However,  in  the  area  of  research,  Dingwall  and 
Scrivener  (1954)  studied  the  stress  distribution  beneath 
slopes  by  the  finite  difference  form  coupled  with  relaxation 
procedures.  They  used  the  theory  of  elasticity  to  determine 
the  shear  and  normal  stresses  for  an  uniform  and  a  rigid 
boundary  embankment. 

Clough  (1960  a,b)  broadened  the  matrix  method  of 
structural  analysis  into  the  general  finite  element  approach 
which  can  be  applied  to  any  structural  mechanics  from  any 
field.  The  matrix  method  of  structural  analysis  was  later 
called  the  finite  element  method. 

Peck  (1967)  pointed  out  that  a  definitive  answer  to  the 
problem  of  progressive  failure  would  require  a  finite 
element  solution  for  a  work-softening  material. 

Duncan  and  Dunlop  (1969)  and  Dunlop  and  Duncan  (1970), 
studied  the  distribution  of  stress  in  and  beneath  slopes  by 
the  finite  element  method.  Their  analyses  were  to  determine 
differences  in  behaviour  of  slopes  in  materials  with  low  and 
high  initial  horizontal  stresses,  representative  of  normally 
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consolidated  and  heavily  over-consolidated  clay  deposits. 

The  material  properties  of  the  slope  were  represented  by 
homogeneous,  linear  elastic,  isotropic  materials. 

On  the  other  hand,  several  analytical  models  were 
developed  to  evaluate  the  mechanistic  approach.  For  example, 
Christian  and  Whitman  (1969)  approach  the  problem  of 
progressive  failure  by  developing  the  differential  equation 
for  displacement  along  the  band  from  equilibrium  of  an 
infinitesimal  element.  Palmer  and  Rice  (1973)  considered  the 
simple  slope  problem  as  an  in-plane  shear  fracture. 

Gibson  (1974)  in  the  Fourteenth  Rankine  Lecture  stated 
that  analytical  methods  draw  attention  to  broad  trends  and 
help  to  distinguish  between  those  factors  that  are  of 
primary  significance  and  those  that  are  of  secondary 
importance . 

Simmons  (1981)  studied  the  behaviour  of  shear  zones. 

The  analyses  of  shearband  yielding  which  involve 
non-weakening  stress-strain  behaviour  are  applied  to  two 
case  histories.  Therefore,  the  future  study  of  the 
analytical  method  should  handle  the  strain-softening 
behaviour  asociated  with  soils  which  are  vulnerable  to 
progressive  failure. 

The  present  research  studies  the  shear  zone  (or  slip 
surface)  behaviour  by  considering  the  strain-weakening 
stress-strain  behaviour.  The  technique  is  applied  to  the 
Edgerton  Slide. 


2.  Description  Of  The  Site  And  Site  Investigation 


2 . 1  General 

The  Slides  occurred  about  48  kilometers  northeast  of 
Wainwright,  Alberta.  The  three  landslides  are  located 
adjacent  to  one  another  and  are  named  the  '  Edgerton  Slides 
'.  The  first  slide  was  termed  the  Edgerton-74  North  Slide. 
The  second  and  third  slides  were  named  the  Edgerton-74  South 
Slide  and  Edgerton-80  Slide,  respectively.  A  plan  view  of 
the  Edgerton  Slides  is  shown  in  Figure  2.1. 

Since  the  first  major  movement  took  place  in  late 
August,  1974,  Thomson  and  Bruce  (1974),  Tweedie  (1976)  and 
Mokracki  (1982)  have  studied  various  aspects  of  the  Edgerton 
Slides . 

Thomson  and  Bruce  (1974)  reported  the  major  features  of 
the  slide  from  a  field  reconnaissance.  Tweedie  (1976)  did  an 
extensive  site  investigation  and  laboratory  testing  program. 
Thomson  and  Tweedie  (1978)  published  a  summary  of  Tweedie' s 
(1976)  work.  Mokracki  (1982)  summarized  the  survey  data  from 
1975  to  1981 . 

2.2  Geology 

Most  of  the  findings  in  this  section  were  obtained  from 
Tweedie  (1976)  and  Mokracki  (1982).  As  Morgenstern  (1977) 
pointed  out  since  the  geological  conditions  in  heavily 
overconsolidated  materials  control  failure  geometry,  some  of 
the  geological  processes  and  materials  will  be  re-emphasized 
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in  order  to  demonstrate  a  part  of  the  reasoning  behind  the 
development  of  the  analytical  procedures  in  Chapter  3. 

2.2.1  Geologic  History 

During  the  period  from  Upper  Cretaceous  to  Lower 
Tertiary,  the  bedrock  formations  were  deposited  in  a 
subsiding  basin  in  Central  Alberta.  Vertical  variation  from 
marine  shale  at  the  base  of  the  Upper  Cretaceous  to 
continental  sandstone  at  the  top  was  common  (  Williams  and 
Burke ,  1964  )  . 

During  late  Mesozoic  and  early  Tertiary  time,  the 
Columbian  and  Pacific  orogeny  transformed  the  Alberta  basin 
from  an  area  of  subsidence  and  deposition  to  one  of  uplift 
and  erosion.  Rutherford  (1928)  estimated  that  about  600 
meters  of  strata  have  been  removed  from  the  study  area 
during  Tertiary  time.  Large-scale  downwasting  and  stagnation 
of  the  Keewatin  ice-sheet,  which  advanced  over  the  area 
during  Pleistocene  time,  modified  the  late  Tertiary 
landscape.  Retreat  of  the  Pleistocene  glaciers  about  10,000 
years  ago,  lead  to  some  topographic  change  due  to  the 
increase  in  river  velocities  and  flow  volumes.  Landslide 
activity  was  started  due  to  the  steep  walled,  post-glacial 
valleys  left  by  the  rapid  downcutting  of  rivers. 

2.2.2  Surficial  And  Bedrock  Geology 


The  glacial  deposits  of  till  are  highly  oxidized  and 
columnar  jointed.  The  average  depth  of  till  within  the  study 


area  is  about  5  meters.  The  composition  of  the  till  is  :  50 
percent  sand,  30  percent  silt,  and  20  percent  clay  sizes  ( 
Bayrock ,  1 967  ) . 
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The  bedrock  of  the  study  area  consists  of  interbedded 
sandstone,  siltstone  and  shale,  and  thin  coal  seams  of  late 
Cretaceous  age  (  Warren  and  Hume,  1939  ).  The  rocks  are 
bentonitic  and  have  a  regional  dip  of  a  meter  per  kilometre 
to  the  southwest. 

The  area  is  underlain  by  the  Bearpaw,  Belly  River  and 
Lea  Park  Formations.  The  Bearpaw  Formation  has  been 
completely  eroded  at  the  study  site.  Thus  the  landslide 
movements  have  occurred  within  the  Belly  River  Formation. 

Possible  ice  shoving  processes  have  occurred  and  were 
observed  over  the  full  height  of  the  scarp  face.  This 
contributed  to  the  brecciated  nature  of  the  bedrock. 

2.3  Description  Of  The  Landslides 

2.3.1  Observations 

The  first  slide  was  discussed  in  detail  by  Thomson 
(1974)  and  Tweedie  (1976).  The  third  slide  was  reported  in 
detail  by  Mokracki  (1982).  The  second  slide  will  be 
discussed  fully  in  this  section,  and  it  will  be  used  for 
later  numerical  analysis. 

Airphotos  of  the  study  area  show  slump  topography  along 
both  sides  of  the  river  valley,  which  indicates  ancient 
landsliding.  Groundwater  discharge  areas  were  found  half  way 
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between  the  river  and  the  local  plain  level.  Toe  erosion  was 
noticed  along  the  river  valley. 

In  the  fall  of  1974,  the  scarp  of  the  Edgerton-74  South 
Slide  was  between  0.45  and  0.6  meters  in  height.  The  slide 
profile  of  1974  will  be  treated  as  pre-slide  profile.  In  the 
summer  of  1982,  the  scarp  on  the  right  flank  was  between  1.2 
and  1.4  meters  in  height.  At  the  same  time,  the  toe  of  the 
south  slide  appears  to  have  cropped  out  at  the  approximate 
location  as  predicted  in  1974. 

2.3.2  Climate 

According  to  records  of  the  Canada  Department  of  Mines 
and  Technical  Services,  1957,  the  climate  of  the  area  is 
sub-humid  continental.  From  the  19  years  of  continuous 
records,  the  average  annual  total  precipitation  is  39.75 
centimeters  of  which  30.84  centimeters  is  rainfall  and  the 
rest  is  89.15  centimeters  snowfall.  Unfortunately,  the 
information  of  the  rainfall  period  (either  long  or  short) 
which  is  of  more  concern  is  not  available. 

2.4  Site  Investigation 

2.4.1  General 

A  review  of  the  past  site  investigation  is  useful  for 
evaluating  the  available  information.  Only  part  of  the 
available  information  can  be  utilized  as  input  data  for 
numerical  analysis.  The  procedural  use  of  the  data  can 


affect  the  numerical  model.  This  will  be  discussed  in  detail 
in  Chapter  3. 

2.4.2  Survey 

The  characteristic  surface  movements  of  the  second 
slide  have  been  monitored  since  early  spring,  1975.  The 
location  of  the  profile  is  shown  in  the  general  plan  of  the 
landslides  (Figure  2.2). 

The  location  of  these  surface  stakes  were  originally 
determined  by  the  changes  in  the  slope  of  the  displaced  mass 
or  by  local  physical  features  along  the  survey  line 

The  typical  recurring  topographic  survey  consists  of 
determining  the  horizontal  distances  and  vertical  elevations 
of  prescribed  control  points  relative  to  a  local  datum.  The 
datum  point  has  been  arbitrarily  selected  and  assigned  an 
elevation  of  300  meters  (local  elevation).  All  horizontal 
distances  to  the  control  point  are  measured  from  this  datum 
point . 

2.4.3  Field  Work  And  Laboratory's  Results 

Tweedie's  (1976)  field  work  consisted  of  subsurface 
exploration  and  instrumentation  programs.  Subsurface 
exploration  included  four  boreholes  and  six  toe  trenches. 
Typical  subsurface  stratigraphy  of  the  second  slide  is  shown 
in  Figure  2.3.  Instrumentation  programs  consisted  of  four 
slope  indicators  and  three  piezometers.  After  a  one  year 
monitoring  program,  all  the  measuring  devices  had  ceased  to 
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function  due  to  continued  slide  movement. 

Laboratory  work  consisted  of  direct  shear  tests. 
Detailed  descriptions  of  the  preparation  of  samples  and  the 
laboratory  program  are  found  in  Tweedie  (1976).  The  shear 
strength  parameters  which  obtained  from  the  laboratory 
program,  are  shown  in  Table  2.1. 

2.5  Section  Summary 

The  preceding  discussion  concluded  that  the  recent 
slides  are,  in  part,  the  reactivation  of  old  landslides. 
Only  surface  movements  have  been  monitored  continuously 
since  May,  1975. 
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pj^gyj-0  2.1  Plan  View  Of  The  Study  Area  With  Slide  Locations 
(Modified  After  Mokracki,  1982) 
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Figure  2.2  General  Plan  Of  Second  Landslide  Showing  Location 
Of  Profiles  And  Boreholes  (Modified  After  Mokracki,  1982) 
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3.  FORMULATION  OF  ANALYTICAL  PROCEDURE 


3 . 1  Introduction 

A  major  portion  of  this  research  work  involved  the 
understanding  and  the  proper  usage  of  the  general  purpose 
linear  and  non-linear  finite  element  analysis  computer 
program  ADINA  for  simulating  the  failure  mechanism  of  the 
landslide . 

This  chapter  can  be  separated  into  two  parts.  The  first 
part  focuses  on  the  exploration  of  the  application  of  ADINA 
in  the  landslide  problems.  Some  of  the  available  options  of 
ADINA  will  be  utilized  for  specific  purposes.  For  example, 
the  option  of  element  death  will  be  used  for  producing  a 
stress-free  boundary  after  excavaton  (or  erosion).  The  three 
material  models  investigated  are: 

a.  isotropic  linear  elastic 

b.  curve  description  with  cracking 

c.  elastic-plastic  (von  Mises,  isotropic  hardening). 
However,  the  analyses  rely  heavily  on  the  isotropic  linear 
elastic  model.  An  attempt  is  made  to  model  stain-softening 
material  behaviour  by  utilizing  the  available  options  of  the 
ADINA  program.  Ultimately,  a  load-transfer  technique  is 
developed  to  model  the  strain-softening  behaviour. 

The  second  part  describes  the  set-up  of  the  numerical 
model  for  finite  element  analysis.  A  simple  model  is  used  to 
study  the  boundary  effect  and  the  appearance  of  the 
softening  zone.  Finally,  an  analytical  procedure  is 
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suggested  for  analysing  the  Edgerton  landslide. 


3.2  Application  Of  Finite  Element  Method  (Program  ADINA)  In 
Modelling  Excavation  And  Material  Softening  Behaviour 

A  three  element  model  will  be  used  to  indicate  the 
applicability  of  the  various  options  of  the  ADINA  program  to 
the  study  of  a  natural  slope.  The  available  options  of 
gravity  loading,  element  death  and  material  models  are 
considered  in  the  following  sections. 


3.2.1  In-Situ  Stresses 

Stresses  within  gravity  structures  due  to  body  forces 
are  of  importance  and  cannot  be  neglected.  According  to  the 
ADINA  user  manual  (Bathe  1976,  section  2.3,  4,  6.1  and 
32.4),  the  mass  proportional  loading  vector  is  established 
from  the  density  of  the  material  and  the  concentrated  mass 
input.  However,  the  aim  of  this  thesis  is  the  study  of  a 
natural  slope;  there  is  no  concentrated  mass  which  derives 
from  a  man-made  structure  to  be  considered  in  the  analysis. 
Hence,  the  gravity  loading  is  directly  proportional  to  the 
density  of  the  material. 

Some  naturally  occurring  sediments  are  deposited  in 
horizontal  layers  where  no  lateral  yielding  occurs.  The 
ratio  of  lateral  to  vertical  stresses  is  known  as  the 
coefficient  of  earth  pressure  at  rest  (Ko).  For  elastic 
isotropic  material,  under  first  loading,  the  value  of  Ko  can 
be  expressed  directly  in  terms  of  Poisson's  ratio  (m);  for 
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example  under  plane  strain  condition, 

u 

Ko  =  - 

1  -  u 

The  in-situ  stresses  are  derived  from  a  switch  on  gravity 
approach;  where  the  vertical  stress  is  due  to  overburden 
load,  and  the  lateral  stress  is  equal  to  Ko  times  the 
vertical  stress. 

Dysli  and  Fontana  (1982)  used  the  ADINA  program  and 
stated  that  the  initial  stress  state  was  created  by  the 
progressive  application  of  gravity  in  about  ten  solution 
steps  in  one  case  and  twenty-two  solution  steps  in  the  other 
case.  However,  the  ADINA  program  is  very  expensive  to  run, 
thus  if  one  solution  step  yields  the  same  results  as  that 
from  many  steps,  instant  gravity  loading  can  be  imposed  on 
the  structure.  For  any  linearly  elastic  analysis,  instant 
gravity  loading  can  be  applied  without  creating  any 
discrepancy  in  finite  element  results.  Therefore,  linearly 
elastic  analysis  is  favoured  for  the  reasons  of  cost  and 
checking  by  hand  calculation  for  a  simple  problem. 

3.2.2  Simulate  Excavation  Using  ADINA 

The  process  of  erosion  is  similar  to  the  process  of 
excavation.  The  differences  are  the  time  scale  and  the 
boundary  conditions.  Erosion  may  take  hundreds  of  years; 
while  excavation  will  take  only  months  or  a  year.  The 
of  excavation  is  defined  by  the  designer;  the  area  of 
erosion  will  be  of  wide  extent. 


area 
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The  element  death  option  is  used  to  simulate  the 
process  of  excavation  (or  erosion).  The  excavated  surface  is 
considered  to  be  stress-free  (Desai  and  Chr i st ian ,  1  977 ) .  The 
stress-free  surface  can  be  created  by  applying  a  set  of 
equivalent  forces  at  nodes  on  the  surface  in  the  direction 
opposite  to  the  direction  of  stresses  due  to  initial  and 
subsequent  loading  conditions. 

The  option  of  element  death  will  yield  a  stress-free 
boundary  under  the  condition  of  no  gravity  loading.  This  is 
shown  by  a  simple  test  which  is  illustrated  in  Figure  3.1. 
However,  under  gravity  loading,  the  program  yields  a  false 
stress-free  boundary  if  the  element  death  option  is  used. 
Kripakov  (1983)  realized  that  the  stiffness  of  the  element 
is  eliminated  at  each  time  step  specified,  but  that  a 
portion  of  the  weight  of  the  element  is  not  effected  (i.e., 
not  eliminated)  if  the  gravity  loading  option  is  used  to 
load  the  structure.  However,  this  is  critical  to  the 
analysis  of  a  natural  slope  as  gravity  is  the  only  loading 
mechanism  which  is  imposed  on  the  slope.  Kripakov  (1983) 
suggested  the  used  of  a  reduced  stiffness  approach  rather 
than  the  death  option  to  simulate  an  excavation  sequence. 

Instead  of  modifying  any  portion  of  the  ADINA  program, 
or  generating  any  complexity  of  the  analysis;  the  problem 
can  be  resolved  by  using  a  thin  layer  of  elements  which  is 
generated  along  the  boundary  of  the  excavation.  Element  1  on 
the  upper  left  of  figure  3.2  can  be  divided  into  two 
elements,  which  leads  to  a  three  element  model.  This  is 
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shown  on  the  upper  right  of  Figure  3.2.  The  reasoning  is 
shown  in  Figure  3.2.  The  error  associated  with  this  approach 
can  be  reduced  by  minimizing  the  thickness  of  the  thin  layer 
(element  1b  in  Figure  3.2).  A  numerical  illustration  of  the 
magnitude  of  the  error  is  shown  in  Table  3.1.  Thus,  using  a 
thin  layer  element  approach,  a  stress-free  boundary  can  be 
obtained . 

NOTE:  the  error,  which  is  generated  from  using  the 
death  option  under  gravity  loading,  does  not  imply  that  the 
program  itself  is  wrong.  This  is  the  standard  finite  element 
procedure  of  distributing  the  weight  of  an  element  to  the 
surrounding  nodes. 

3.2.3  Material  Models 

Of  the  fourteen  material  models  available  in  ADINA, 
only  three  are  considered  in  this  research.  These  are: 

a.  curve  description  with  cracking 

b.  elastic-plastic  (von  Mises,  isotropic  hardening) 

c.  isotropic  linear  elastic. 

Considering  the  curve  description  model,  in  the 
beginning  some  ADINA  users  thought  that  the  curve 
description  model  can  be  used  for  modeling  the  post-failure 
behaviour  of  strain  softening  material.  A  decrease  of 
strength  occurs  after  shear  straining  past  the  peak  value, 
which  causes  a  progressive  type  of  failure  in  stiff, 
fissured  clay  slopes.  The  stress-strain  curve  of  a  general 
strain-softening  material  is  shown  in  Figure  3.3.  If  one 
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uses  the  curve  description  model  for  modeling  the 
post-failure  behaviour  of  a  strain-softening  material,  a 
matric  inversion  difficulty  in  ADINA  is  created.  Both 
Kripakov  (1983)  and  Dysli  (1982)  have  already  expressed 
their  scepticism  concerning  this  function  of  ADINA. 

If  the  stresses  exceed  the  yield  stress,  the  results  of 
an  elastic-plastic  analysis  will  indicate  the  location  of 
the  plastic  zone.  If  the  stresses  are  less  than  the  yield 
stress,  the  results  of  such  an  analysis  are  identical  to 
those  from  a  linear  elastic  analysis.  The  elastic-plastic 
model  cannot  be  used  for  a  strain-softening  material. 

In  usual  engineering  analysis,  a  non-linear  analysis  of 
a  program  is  always  preceded  by  a  linear  analysis.  The 
advantages  of  an  isotropic  linear  elastic  analysis  are: 

1.  the  results  can  easily  be  checked  by  hand  calculations 
for  some  cases. 

2.  the  least  number  of  input  parameters  are  required  for 
the  analysis. 

3.  the  strain-softening  material  behaviour  cannot  be 
modelled  by  any  available  material  models.  Therefore, 
there  is  no  real  advantage  to  using  any  sophisticated 
material  model  at  this  time. 

Hence,  material  models  of  both  elastic-plastic  and  curve 
description  with  cracking  are  not  considered  in  the 
following  analyses. 

At  first  glance,  material  softening  can  be  modelled  by 
using  both  the  options  of  element  birth  and  element  death 
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simultaneously  within  the  ADINA  program.  The  approach  is 
similar  to  the  incremental  elasticity  approach  in  which  the 
Young's  Modulus  (E)  is  decreasing.  The  stress-strain  curve 
of  the  incremental  elasticity  approach  is  shown  in  Figure 
3.4.  The  sequence  of  operation  is  that  element  1b  (refer  to 
Figure  3.2)  with  El  is  killed  and  is  replaced  by  a  new 
element  having  E  corresponding  to  E2  (refer  to  Figure  3.4). 
Process  of  replacing  element  1b  with  an  element  of  changing 
E  according  to  Figure  3.4  is  continued.  As  strain  increases, 
E  decreases.  Hence,  strain-softening  can  be  modelled.  On  a 
theoretical  basis,  this  approach  is  better  for  modelling  the 
strain-softening  behaviour  than  any  available  material 
models.  However,  if  the  death  option  is  used  in  the  first 
step,  the  stresses  within  the  element  will  turn  to  zero.  If 
the  birth  option  is  used  in  the  second  step,  the  element 
will  carry  no  initial  stresses.  Therefore,  the  ADINA  can  be 
used  to  perform  the  incremental  elasticity  approach,  but  the 
stress  output  will  be  zero  for  the  second  step  of  the 
operation  of  death  and  birth  sequence  of  the  program  used. 

The  previous  discussions  illustrate  that  the  ADINA 
program  should  be  used  in  conjunction  with  another  program 
in  order  to  model  the  strain-softening  behaviour.  Of  course, 
the  best  approach  is  to  modify  the  ADINA  program  so  that 
this  function  is  built  into  the  ADINA  program  itself.  The 
technique  of  the  new  approach  will  be  discussed  thoroughly 
in  the  next  section,  but  the  work  of  modifying  the  ADINA 
program  to  include  the  new  material  model  is  left  for  future 
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research. 

3.3  Load-Transfer  technique 

One  of  the  necessary  preconditions  for  progressive 
failure  to  occur  is  that  the  material  of  a  slope  must 
display  strain-softening  behaviour.  The  load-transfer 
technique,  to  be  described  in  the  following  paragraphs,  is 
actually  a  strain-weakening  approach.  While  this  is  not  the 
same  as  strain-softening,  it  is  a  better  approach  than  most 
others  for  understanding  the  progressive  failure  process. 

3.3.1  Background  Information 

The  technique  is  analogous  to  the  stress  transfer 
method  (Zienkiewicz  and  et.al.,  1968).  The  latter  method  is 
used  for  the  stress  analysis  of  a  rock  mass  which  cannot 
sustain  tension.  The  load-t ransf er  technique  is  used  for 
reducing  shear  modulus  of  a  shear  zone  material  which  cannot 
sustain  excessive  shear  stresses.  The  assumed  load 
redistribution  approach  is  not  capable  of  reproducing  the 
true  strain-softening  behaviour,  however,  it  does  provide 
useful  method  for  understanding  the  stability  of  a  natural 
slope . 

3.3.2  Description 

This  technique  utilizes  the  available  program  ADINA 
within  the  Department  of  Civil  Engineering  at  the  University 
of  Alberta  and  the  load  redistribution  program. 
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The  essential  steps  of  this  approach  can  be  described 

as  follows: 

1.  analyse  the  problem  as  an  isotropic  elastic  case  using 
the  ADINA  program.  The  loading  mechanism  is 
gravitational  force. 

2.  reduce  the  elastic  modulus  at  certain  locations  (eg. 
shear  zone  or  slip  surface)  which  may  not  be  capable  of 
sustaining  excessive  shear  stress. 

3.  the  restraining  forces  are  generated  due  to  the 
reduction  of  shear  strength  in  terms  of  elastic  modulus. 
These  forces  are  obtained  from  the  load-transfer 
program.  Total  stress  (ot)  of  the  element  can  be 
separated  into  two  components. 


am  is  the  stress  that  can  be  sustained  by  the  element 
with  the  reduced  modulus. 

o n  is  the  excess  stress  that  cannot  be  sustained  by  the 
element  due  to  softening  which  must  be  redistributed  to 
other  parts  of  the  structure. 

Therefore,  the  new  stress  for  the  element  will  be 
am  and  the  equivalent  load  (R)  from  that  must  be  applied 
to  redistribute  this  excess  stress. 


More  detail  discussion  of  this  procedure  is  given  in 
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Appendix  A. 

4.  the  application  of  forces  to  relieve  the  restraining 
forces  is  used  to  maintain  equilibrium.  The  ADINA 
program  is  re-utilized  again.  The  results  will  be  the 
incremental  stresses. 

5.  the  stresses  are  computed  in  such  a  way  that  the 
incremental  stress  (6a)  due  to  the  applied  load  (R)  will 
be  superimposed  on  am,  but  not  on  at. 

The  flow  chart  is  shown  in  Figure  3.5. 

3.4  Examples  Of  The  Load-Transfer  Program 

Two  examples,  namely  pure  shear  and  bending  of  a  beam 
by  uniform  load,  are  used  to  verify  the  load-transf er 
program  and  to  demonstrate  the  technique  of  reducing  shear 
modulus.  The  data  input  instruction  and  the  source  code  of 
the  program  is  shown  in  the  Appendix  B. 

3.4.1  Pure  Shear 

The  pure  shear  model  and  its  finite  element 
idealization,  the  material  properties  and  the  loading 
conditions  are  shown  in  Figure  3.6.  If  the  shear  modulus  (G) 
is  reduced  from  385  kPa  to  96  kPa,  the  generated  loads  due 
to  excess  shear  stress  have  to  be  redistributed  to  the  rest 
of  the  structure.  Theoretically,  the  stresses  and  the 
strains  from  the  analysis  with  the  material  property  of  96 
kPa  should  be  the  same  as  the  summation  of  the  stresses  from 
that  of  385  kPa  and  from  excess  shear  stress. 
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Refer  to  Figure  3.6, 

^xyk  —  ^  x  y  i  5  C  x  y 

7k  =  7 i  +  ^7 

where 

axyk  and  7k  =  final  shear  stress  and  shear 
strain  respectively 

6axy  and  67  =  calculated  incremental  shear 
stress  and  shear  strain  respectively 
axyi  and  7  j  =  initial  shear  stress  and  shear 
strain  respectively 
For  the  case  1 : 

-in  step  1,  the  stresses  and  strains  are  calculated 
from  the  ADINA  program  with  the  Young's  Modulus  of 
1000  kPa . 

—in  step  2,  the  excess  shear  stresses  in  terms  of 
loads  are  calculated  from  the  load-transfer  program 
with  the  reduction  of  the  Young's  Modulus  from  1000 
kPa  to  250  kPa. 

—in  step  3,  the  excess  loads  from  step  2  are  applied 
to  the  ADINA  program  where  the  incremental  stresses 
and  strains  are  calculated  and  are  superimposed  on 
those  in  step  2. 

If  the  shear  modulus  is  reduced  to  a  quarter  of 
its  original  magnitude,  the  strain  has  to  increase 
in  order  to  reach  the  same  stress  level.  Figure  3.7 
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shows  the  stress  path  of  the  redistribution  of 
stress  due  to  the  reduction  of  shear  modulus.  If  the 
material  cannot  take  the  original  stress  for 
whatever  reasons,  the  material  properties  (in  terms 
of  Young's  Modulus  or  Poisson's  Ratio)  have  to  be 
reduced  to  a  lower  value.  The  excess  shear  stress 
will  be  developed  due  to  the  reduction  of  Young's 
Modulus.  Therefore,  the  redistribution  of  the  excess 
shear  stress  will  generate  an  additional  strain. 

If  the  model  is  analysed  by  using  the  ADINA 
program  with  the  Young's  Modulus  of  250  kPa ,  (i.e. 
case  2,  refer  to  Figure  3.6)  the  results  (both  axyk 
and  7k )  of  case  2  should  be  the  same  as  those  of 
case  3  in  step  1.  Therefore,  it  is  concluded  that 
the  load-t ransf er  program  is  acceptable. 

3.4.2  Bending  Of  A  Beam  By  Uniform  Load 

This  example  serves  two  purposes.  These  are: 

a.  apply  the  load-transfer  program  to  a  complex 
problem. 

b.  compare  the  result  from  the  new  technique  with  the 
closed  form  solution  from  the  theoretical  approach. 

A  cantilever  beam  subjected  to  uniform  distributed  load 
is  being  considered.  The  finite  element  idealization  of  the 
problem  and  the  material  properties  are  shown  in  Figure  3.8. 
The  material  property  (primarily  elastic  modulus)  is 
reduced,  and  the  Poisson's  ratio  is  kept  constant. 
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The  calculation  of  load-redistribution  is  based  on  the 
plane  strain  condition;  however  the  theoretical  solution  for 
this  problem  is  based  on  simple  beam  theory  (plane  stress). 
The  equation  which  is  used  for  calculating  the  deflection  at 
node  35  (refer  to  Figure  3.8)  is: 

w  1 4 

6  =  - 

8EI 

where 

w  =  uniformly  distributed  load 
1  =  length  of  the  beam 
E  =  Young  Modulus 
I  =  moment  of  the  inerita 
6  =  deflection 

The  result  from  the  load  redistribution  approach 
follows  closely  the  one  from  the  theoretical  solution.  The 
results  are  presented  in  Figure  3.9. 

3.5  Summary  Of  The  Load-Transfer  Technique 

The  new  technique  has  been  described  and  proven  to  be 
in  the  working  stage.  The  technique  is  actually  an 
additional  material  model  for  analysing  any  strain  weakening 
problems.  The  last  two  examples  demonstrate  how  the  shear 
modulus  was  reduced  for  the  entire  structure. 

The  technique  of  reduction  of  the  shear  modulus  for  a 
part  of  the  structure  will  depict  the  picture  of  * 
non-homogeneous  behaviour.  Additionally,  for  a  stability 
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problem,  the  shear  strength  will  be  mobilized  along  the  slip 
surface  or  a  pre-existing  weak  zone  (e.g.  shear  zone).  The 
next  section  will  study  this  question  in  detail  and  propose 
a  method  for  analysing  the  Edgerton  Slide. 

3.6  Behaviour  Of  Simple  Model 

Usually,  erosion  of  a  landscape  proceeds  in  a 
systematic  way  so  that  landforms  evolve  through  a  series  of 
stages  in  which  the  ultimate  landscape  is  reduced  to  a 
surface  of  low  relief  (Hamblin  and  Howard,  1975). 

The  simple  model  which  is  shown  in  Figure  3.10 
undergoes  the  following  steps: 

1.  assign  material  properties  to  the  whole  structure 

2.  switch-on  gravity 

3.  excavate  part  of  the  structure 

4.  reduce  the  modulus  in  shear  zone. 

The  ADINA  approach  requires  a  pre-existing  shear  zone 
before  the  excavation  process;  while  the  present 
load-transfer  approach  allows  the  softening  process 
after  the  excavation  stage,  (refer  to  Figure  3.10) 

5.  re-analyse  the  problem  with  new  modulus  using  ADINA. 

The  aim  of  the  simple  model  is  to  determine  how  much 
variation  of  the  result  will  be  arised  from  the  modelling 
techniques . 
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3.6.1  Introduction 

Theoretically,  progressive  failure  indicates  the 
spreading  of  the  failure  over  the  potential  surface  of 
sliding  from  a  point  or  a  line  toward  the  boundaries  of  the 
surface  or  vice  versa.  Therefore,  the  strength  properties 
across  the  shear  zone  will  not  be  uniform.  Two  questions 
should  be  answered  in  order  to  understand  the  failure 
mechanism.  These  questions  are: 

1.  how  much  reduction  of  shear  modulus  (in  term  of  Young's 
Modulus)  is  required  to  induce  the  stress  concentration? 

2.  where  will  the  softening  zone  be  first  initiated  and  in 
what  direction  does  this  softening  zone  progress? 

Several  people  have  been  investigating  the  second 
question  for  a  long  time,  for  example,  Palmer  and  Rice 
(1973)  and  Chowdhury  (1978). 

The  research  of  this  thesis  involves  only  the  first 
question.  Although  the  second  question  is  studied,  attempts 
have  been  made  with  no  conclusive  result  which  is  required 
on  fracture  mechanics  (J-integral)  can  be  drawn  from  this 
study.  The  variation  of  the  results  should  arise  from  the 
change  of  shear  strength  within  the  shear  zone.  Hence,  it  is 
important  that  the  boundary  effect  not  influence  the  result. 

3.6.2  Boundary  Effect 

The  aim  of  this  section  is  to  determine  how  much 
variation  will  be  derived  from  changing  the  boundary 


conditions . 
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The  nodes  at  A  and  B  are  assumed  to  be  fixed,  (refer  to 
Figure  3.10)  This  is  due  to  the  stability  of  the  structure. 
The  nodes  along  the  upslope  and  downslope  side  can  be 
assumed  to  be  on  rollers  because  deformation  will  take  place 
in  a  vertical  direction  under  gravity  loading.  However,  the 
nodes  along  the  bottom  boundary  can  be  assumed  to  be  either 
fixed  or  on  rollers.  If  the  boundary  is  set  far  enough  away, 
the  results  from  either  fixed  or  roller  bottom  boundary 
should  be  approximately  same. 

It  is  assumed  that  the  whole  structure  is  uniform  at 
this  moment.  The  material  properties  are: 

1)  Young’s  Modulus  =  137,900  kPa 

2)  Poisson  Ratio  =  0.42 

The  sequence  of  the  analysis  consists  of  switch-on 
gravity  and  then  the  excavation  process.  Therefore,  the  only 
variable  for  this  problem  comes  from  the  geometry. 

The  accuracy  can  be  increased  by  increasing  the  number 
of  nodes  and  elements.  However,  there  is  a  limit  to 
increasing  these  two  parameters  because  the  computation  time 
follows  a  power  rule  of  the  number  of  nodes  and  elements.  An 
number  of  trials  involving  various  meshes  was  done  to 
accommodate  the  features  of  a  stress-free  boundary  and  shear 
zone.  The  final  configuration  of  the  mesh  used  in  the 
analyses,  consists  of  280  nodes  and  255  4-node  elements 
(refer  to  Figure  3.11(b)  ). 

The  results  derived  from  roller  bottom  boundary  differ 
from  those  derived  from  fixed  bottom  boundary.  This 
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difference  reduces  as  the  dimensions  of  B  and  H  (refer  to 
Figure  3.10)  increase.  However,  the  difference  can  be 
narrowed  down  to  10  to  20  percent  without  substantially 
increasing  the  dimension  of  B  and  H.  Figure  3.11  shows  two 
of  the  possibilities.  Table  3.2  shows  the  percentage 
difference  of  the  displacements  obtained  from  both  the  fixed 
and  roller  bottom  boundary. 

When  a  relatively  hard  material  is  absent,  it  is 
necessary  to  establish  finite  boundaries  within  which  the 
results  will  not  be  changed  due  to  the  boundary  effect.  It 
is  assumed  that  the  boundary  effect  of  the  configuration  of 
Figure  3.11b  is  minimal. 

3.6.3  Softening  Zone 

The  schematic  illustration  of  the  simple  model  is  shown 
in  Figure  3.12.  The  simple  model  consists  of  four  different 
layers.  Their  properties  are  shown  in  Table  3.3.  The 
thinnest  layer  is  assumed  to  be  the  shear  zone  (or  slip 
surface).  As  the  material  within  the  shear  zone  will  undergo 
non-uniform  strength  weakening,  certain  locations  within  the 
shear  zone  will  have  lower  strength  properties.  Therefore, 
artifical  reduction  of  Young’s  Modulus  to  either  one  or  five 
percent  of  its  original  value  is  assigned  for  either 
location  A  or  B.  (refer  to  figure  3.12) 

The  purpose  of  this  analysis  is  to  determine  how  the 
stress  variation  will  take  place  if  the  strength  within  the 
shear  zone  is  reduced. 
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3.6.3.  1  Procedures 

The  procedures  are  shown  in  Figure  3.10;  namely, 
ADINA  approach  and  load-transfer  approach.  However,  the 
concepts  of  the  two  approaches  are  different.  The  former 
approach  assumes  that  the  softening  zone  exists  before 
the  excavation  process.  The  latter  approach  assumes  that 
the  softening  zone  is  generated  after  the  excavation 
process . 

The  operational  sequence  of  the  simple  numerical 
model  consists  of  two  steps.  The  boundary  of  the  first 
step  is  MNYX  as  shown  in  Figure  3.11b  which  resembles 
the  flat  and  horizontal  landscape.  The  boundary  of  the 
second  step  is  MABCYX  as  shown  in  Figure  3.11b,  which 
portrays  the  eroded  (or  excavated)  landscape.  The  only 
difference  between  these  two  approaches  comes  from  the 
assigning  of  the  material  properties  to  each  of  the 
elements  in  the  first  step.  For  the  ADINA  approach,  the 
material  properties  remain  unchanged.  On  the  other  hand, 
the  load-transfer  approach  assumes  that  the  shear  zone  • 
material  is  uniform  in  the  first  step.  The  material 
properties  within  a  given  area  of  the  shear  zone  will  be 
reduced  after  the  second  step.  The  detailed  description 
has  been  mentioned  in  Section  3.3.2. 

Although  the  final  material  properties  derived  from 
both  approaches  are  same,  the  cost  of  a  run  using  the 
ADINA  approach  is  approximately  half  of  that  using  the 
load-transfer  approach.  Will  the  results  derived  from 
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these  two  approaches  be  significantly  different? 

Both  approaches  are  used  in  analysing  several 
trials.  These  trials  are  shown  in  Table  3.4.  These 
trials  may  indicate  how  sensitive  the  result  will  be  due 
to  a  reduction  of  the  strength  in  terms  of  Young's 
Modulus . 

3. 6. 3. 2  Observation  And  Comparison  Of  Results 

The  output  from  any  finite  element  program  will  be 
in  terms  of  displacement  (strain)  and  stress. 

DISPLACEMENT 

A  typical  pattern  of  the  displacement  arrows  due  to 
the  excavation  process  (or  valley  erosion)  is  shown  on 
Figure  3.13.  The  results  of  the  surface  displacements 
from  the  ADINA  approach  and  the  load-transfer  approach 
are  shown  in  Figure  3.14  for  trial  1.  (refer  to  Table 
3.4)  and  in  Figure  3.15  for  trials  2  of  the 
load-transfer  approach  and  3  of  the  ADINA  approach.  The 
phenomenon  of  valley  rebound  (Matheson,  1973)  can  be 
shown  by  the  surface  displacements. 

From  a  comparison  of  the  results  of  the  analyses, 
the  load-transfer  approach  tends  to  yield  a  lower 
surface  rebound  value  than  that  of  ADINA.  The  lower  the 
Young's  Modulus  (eg.  trial  1  refer  to  Table  3.4)  yields 
a  lower  rebound  value.  Both  the  location  and  the  size  of 
the  softening  zone  will  affect  the  value  of  surface 
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rebound  value.  However,  the  size  of  the  softening  zone 
is  proportional  to  the  value  of  the  surface  rebound. 

STRESSES 


The  only  loading  mechanism  of  the  first  step  in  the 
analysis  is  the  gravitational  force.  Since  the  simple 
model  consists  of  uniform  layers  and  no  horizontal  shear 
stress,  therefore,  the  vertical  stress  is  the  minor 
principal  stress  (sign  convention  is  that  compression  is 
negative).  The  horizontal  stress  should  be  the  major 
principal  stress.  However,  the  analytical  results  show  a 
small  discrepancy  (about  1  percent)  with  the  closed-form 
solutions . 

Four  uniform  layers  as  the  case  of  step  1  of  the 
load-transfer  approach  will  yield  zero  horizontal  shear 
stress  within  the  structure.  The  horizontal,  vertical 
and  maximum  shear  stress  contours  will  be  a  number  of 
horizontal  lines. 

For  the  ADINA  approach,  the  layer  three  (refer  to 
Figure  3.12)  is  not  uniform  in  step  1.  Therefore,  the 
concentration  of  the  shear  stress  contours  is  primarily 
due  to  the  non-homogeneous  effect.  This  effect  appears 
more  if  the  softening  zone  is  located  at  A  (refer  to 
F i gure  3.12). 

In  the  second  step,  a  stress  discontinuity  will  be 
obtained  at  the  softening  zone.  The  overall  pattern  of 
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stress  contours  is  almost  identical  from  either  one  of 
the  approaches.  The  softening  zone  at  location  A  seems 
to  have  a  large  stress  concentration  area.  This  can  be 
illustrated  by  several  shear  stress  (rxy)  contours  as 
shown  on  Figures  3.16  and  3.17. 

However,  if  the  results  are  studied  in  depth,  one 
can  conclude  that  the  load-transfer  approach  yields  a 
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2.  the  load  transfer  approach  yields  uniform  stress 
across  the  softening  zone. 

3.  the  ADINA  approach  yields  a  relative  peak  stress  at 
the  mid  -section  of  the  softening  zone. 

Additionally,  it  seems  that  a  large  reduction 
of  elastic  parameters  is  required  to  induce  a  stress 
concentration.  From  the  experience  of  this  analysis, 
the  Young's  Modulus  has  to  be  reduced  to  one  to  five 
percent  of  its  original  value  to  observe  a 
noticeable  stress  concentration. 

Although  the  results  do  not  show  any 
significant  difference  between  the  two  approaches, 
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it  is  important  to  recognize  that  the  load-transfer 
approach  matches  the  geological  failure  mechanism 
sequence  of  the  Edgerton  Slide.  As  noted  previously, 
the  load-transfer  method  does  not  assume  a 
pre-existing  weak  zone.  Therefore,  the  load-t ransf er 
approach  will  be  used  for  the  analysis  in  Chapter  4. 

3.7  Chapter  Summary 

The  preceding  discussion  summarized  most  of  the 
required  techniques  which  will  be  used  for  the  deformation 
analysis  of  the  Edgerton  Slide.  The  correct  usage  of  the 
computer  program  of  ADINA  was  emphasized  so  that  meaningful 
results  can  be  obtained.  The  prime  purpose  of  this  chapter 
is  to  provide  a  conceptual  feeling  for  the  results  so  that 
fewer  trials  will  be  required  in  Chapter  4  and  hence  the 
computational  cost  will  be  reduced. 
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Table  3.1  Error  Associated  With  The  Thin  Layer  Approach 


Assume 
H  =  10  m 

P  =  2140  kg/cu.  m 


h/H 

VERTICAL 

STRESS  ( kPa ) 

(at  the 

centre  of  element 

0.05 

110.1030 

0.01 

105.9086 

0.005 

105.3843 

These  are  compared  to  the 

Ideal  Case  where 

104.8600 

No  Thin 

Layer 

209.7200 

(  i  .  e  .  2 

element 

model ) 

Note  : 

1.  h/H  cannot  equal  zero. 

2.  h/H  =  .005  used  in  calculation  of  this 
comparison . 

3.  percentage  error  of  h/H  =  .005  is  0.5  percent. 
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Table  3.2  Comparison  Of  Displacement  Output 


CONDITION 

DISPLACEMENT  OUTPUT  FOR  NODES 

99 

108 

* 

FREE 

0.0210 

0.1409 

FIGURE  3.11a 

FIXED 

0.0313 

0.1169 

(33%  higher) 

(17%  lower ) 

FREE 

0.0714 

0.2305 

FIGURE  3.11b 

FIXED 

0.0728 

0.2055 

( 1 9%  higher) 

(  17%  lower) 
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Table  3.4  Location  And  Material  Properties  Of  The  Trials 


TRIAL 

REDUCED 

YOUNG’S  MODULUS 

(MPa) 

(original  lOOMPa) 

LOCATION 

(refer  to  Figure  3.12) 

1 

1 

A 

2 

5 

A 

3 

1 

B 

4 

5 

B 
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STEP  1 


LOAD  THE  STRUCTURE  AS  SHOWN 


STEP  2 


KILL  (DESTROY)  ELEMENT  1 


2 

TTrr  7T77 

CONCLUSION  : 


THE  STRESS  OF  ELEMENT  2  SHOULD 
BE  ZERO  IF  THE  STRESS-FREE  BOUNDARY 
IS  VALID. 


Figure  3 
Gravity 


.1  A  Simple  Test  Of  Element  Death  Option  Without 
Loading  (For  Simplicity,  Only  Two  Elements  Are  Used) 
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ELEMENT 


FROM  FIGURE  3 . 1 


MODIFIED  STRUCTURE 


STEP  1  LOAD  THE  RIGHT  HAND  STRUCTURE  GRAVITY 

STEP  2  KILL  (DESTROY)  ELEMENTS  l.a  AND  l.b 

The  vertical  stress  of  element  2  (a) 
should  be  : 

o  =  p  g  hd  -  ( 1 ) 

However,  the  internal  calculation  of 
ADINA  program  would  be  : 

h 

a=pghd+pg-  -  (2) 

2 

As  h  approaches  zero,  a ( 2)  approaches  o( 1) 

Figure  3.2  A  Simple  Test  Of  Element  Death  Option  With 
Gravity  Loading 
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SHEAR 


STRESS 


VOLUME 

STRAIN 


Figure  3.3  The  Strain-Softening  Model 


STRESS 
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NOTE  :  The  value  of  E  is  the  slope  of 


the  stress-strain  curve. 


Figure  3.4  The  Incremental  Elasticity  Model 
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Figure  3.5  Flow  Chart  Of  The  Load-Transfer  Technique 
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1  0  m 


PURE  SHEAR  MODEL 


FINITE  ELEMENT  IDEALIZATION 


MATERIAL  PROPERTIES  : 

Density  (p)  =  0  kg/cu.m. 

Poisson's  Ratio  in)  =  0.3 
Young's  Modulus  (El)  =  1000.0  kPa 
Young's  Modulus  (E2)  =  250.0  kPa 
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96 

289 
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0.45 
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6.24 
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0.0 
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0.0 

6.24 

6.24 

Figure  3.6  Example  Of  Pure  Shear 


SHEAR  STRESS  (kPa) 
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SHEAR  STRESS  VS  STRAIN 


Figure  3.7  Stress  Path  For  The  Reduction  Of  Shear  Modulus 
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Example  3.8  Example  Of  Bending  Of  A  Beam  By  Uniform  Load 


DEFLECTION  RT  NODE  35  (m) 
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Figure  3.9  Deflection  vs  Young's  Modulus  For  The  Bending  Of 
A  Beam  By  Uniform  Load 
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ups  1  ope 
s  i  de 


downs  lope 
s  i  de 


ADINA  APPROACH  LOAD-TRANSFER  APPROACH 


Figure  3.10  The  Simple  Model  With  1.  ADINA  approach  2. 
Load-Transfer  approach 
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Figure  3.11  Meshes  Of  The  Simple  Model 
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SCHEMATIC  DIAGRAM  OF  THE  SIMPLE  MODEL 


MAGNIFICATION  OF  THE  SHEAR  ZONE  (element  numbering) 


LOCATION  OF  A  INCLUDES  ELEMENTS  179,  180,  181  &  182 

LOCATION  OF  B  INCLUDES  ELEMENTS  184,  185  &  186 


Figure  3.12  Schematic  Illustration  Of  Shear  Zone 


ELEVATION  C  M  ) 
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Figure  3.13  A  Typical  Pattern  Of  The  Displacement  Arrows  Due 
To  Excavation  Process  por  Trial  3  (ADINA  Approach) 
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Figure  3.14  a)  Surface  Displacements  For  Trial  1  (ADINA 
Approach ) 


Figure  3.14  b)  Surface  Displacements  For  Trial  1 
(Load-Transfer  Approach) 
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Figure  3.15  a)  Surface  Displacements  For  Trial  2 
(Load-Transfer  Approach) 


Figure  3.15  b)  Surface  Displacements  For  Trial  3  (ADINA 
Approach ) 


3:  x  .  3  ;  :  3. 1  *ir: 
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Figure  3.16  a)  Shear  Stress  (rxy)  Contours  For  Trial  1 
(ADINA  Approach) 


O 


Figure  3.16  b)  Shear  Stress  (rxy)  Contours  For  Trial  1 
(Load-Transfer  Approach) 
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Figure  3.17  a)  Shear  stress  (rxy)  Contours  For  Trial  3 


(ADINA  Approach) 


O 


Figure  3.17  b)  Shear  Stress  (rxy)  Contours  For  Trial  3 
(Load-Transfer  Approach) 
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Figure  3.18  Illustration  Of  The  Locations  Of  The  Sections 
Used  For  Figures  3.19  and  3.20 


ADINA  Approach 
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Figure 


3.20  The  Change  Of  The  Maximum  Shear  Stress  Acros 


Several  Sections  For  Trial  3 


4.  FINITE  ELEMENT  ANALYSIS  OF  EDGERTON  SLIDE 


4.1  Aim  Of  The  Analysis 

This  chapter  presents  an  example  of  the  application  of 
the  analytical  method  to  a  field  problem.  A  deformation 
analysis  of  the  Edgerton  Slide  was  carried  out  for  this 
research . 

The  deformation  analysis  is  able  to  clarify  the  nature 
of  the  failure  mechanism  as  will  be  shown  in  this  chapter. 
This  has  some  important  implications  for  stability  design 
which  is  usually  based  on  limit-equilibrium  approach. 

The  relative  importance  of  the  following  variables 
should  be  established  before  generalizing  slopes  on  the 
basis  of  the  geometry,  scale  and  soil  properties.  These 
factors  are  quoted  from  Bishop  (1971) 

a.  the  relationship  between  post-peak  drop  off  in 
strength  and  displacement. 

b.  the  swelling  characteristics  of  the  soil. 

c.  the  prepeak  stress  deformation  characteristics 
of  the  soil  under  the  appropriate  conditions  of 
stress  change. 

d.  the  value  of  the  coefficient  of  earth  pressure 
at  rest  before  the  formation  of  the  slope. 

e.  the  geometry  and  scale  of  the  slope. 

f.  the  long  term  flow  pattern  of  the  ground  water. 

All  these  factors  will  be  considered  except  points 

(b)  and  (f);  in  an  attempt  to  correlate  the  field  data 
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with  the  analytical  results.  The  ultimate  goal  is  to 
understand  a  failure  mechanism  and  a  secondary 
deformation  pattern  which  will  be  used  as  a  guide  to 
assess  the  stability  of  a  natural  slope  » 

The  logic  behind  the  analysis  is  of  considerable 
importance  in  extrapolating  from  experience  from  the 
Edgerton  slide  to  other  similar  slides. 

The  dilatancy  which  accompanies  failure  and 
post-failure  behavior  of  over-consolidated  clay,  is 
localized  in  a  very  narrow  zone  around  the  failure 
surface.  This  behavior  is  inadequately  represented  in 
most  of  the  numerical  models  recently  proposed  for 
describing  the  strain-softening  and  dilatancy  of  stiff 
clays . 

The  purpose  of  this  chapter  is  to  gain  an  insight 
into  the  landslide  mechanism.  The  general  purpose  finite 
element  program  coupled  with  the  load-transfer  program 
is  used  to  predict  the  stability  of  st i f f - f i ssured  clay 
slopes.  The  process  of  modelling  the  landslide  from 
beginning  to  end  involves  large  displacement  finite 
element  formulation.  Additionally,  the  problems  of 
rupture  and  cracking  are  very  difficult  to  model.  It 
seems  that  none  of  the  existing  programs  can  handle 
these  problems.  Therefore,  the  work  of  modelling  the 
whole  process  of  any  landslide  is  left  for  future 


research . 
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4.2  Possible  Mode  Of  Failure 

Thomson  and  Tweedie  (1977)  postulated  that  the  failure 
of  the  Edgerton  Slide  occurred  due  to  a  gradual  loss  of  soil 
strength,  manifested  by  a  virtual  disappearance  of  cohesion, 
with  the  final  triggering  mechanism  being  a  spring  time  rise 
in  the  pore  pressure  within  the  slide  mass.  The  existence  of 
the  pre-sheared  failure  plane  is  the  result  of  several 
earlier  stages  of  landslide  activity. 

The  following  analyses  will  rely  heavily  on  the  choice 
of  the  strength  parameters.  Both  the  time  element  and  the 
piezometric  level  are  not  considered  in  the  analysis.  A 
weaker  strength  parameter  is  assigned  to  the  material  below 
the  water  table.  Even  though  delayed  pore  pressure 
equalization  is  an  important  factor  in  an  analysis  of  slope 
instability,  it  is  most  likely  that  the  previous  slide 
history  or  pore  pressures  are  not  known  in  detail. 

4.3  Field  Work  Essential  For  An  Analysis 

The  movement  of  the  Edgerton  Slide  was  carefully 
monitored  by  Tweedie  (1975)  using  three  slope  indicators 
along  the  slide  profile.  The  location  of  these  indicators  is 
shown  in  Figure  2.2.  The  slope  indicator  data  are  shown  in 
Figures  4.1  to  4.3.  Yearly  surface  movements  had  been 
documented  by  Mokracki  (1982). 

The  slide  profile  has  been  monitored  since  1975.  From 
the  interpretation  of  the  yearly  surface  movements,  the  mass 
movement  can  be  divided  into  four  individual  blocks. 
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Therefore,  three  probable  rupture  surfaces  will  likely  be 
presented.  These  rupture  surfaces  will  be  considered  in  the 
design  of  the  finite  element  mesh,  which  is  shown  on  Figure 

4.4. 

From  field  observations,  a  consistent  cracking  pattern 
is  found  throughout  the  site;  particularly  the  area  close  to 
the  bulging  and  the  toe.  The  hexagonal  cracking  pattern  is 
observed  in  the  field  and  is  difficult  to  model  by  finite 
element  analysis.  However,  the  surficial  cracking  pattern 
may  not  lead  to  catastrophic  failure.  This  cracking  pattern 
may  only  indicate  the  area  of  tension. 

4.4  Uncertainty 

This  portion  of  the  work  gained  from  personal 
discussions  with  Dr.  John  Hutchinson  during  the  fall  of 
1982.  As  the  major  work  of  this  thesis  was  to  model  the 
natural  slope  by  the  finite  element  method,  a  few  possible 
conditions  can  occur  as  a  result  of  limited  field 
information.  These  are  : 

1.  the  relative  displacements  between  the  various  layers 
are  uncertain. 

The  problem  can  be  narrowed  down  to  the  shear  band 
problem.  Within  the  shear  zone,  the  velocity  gradient 
can  be  expected  to  be  higher  than  those  of  the  adjacent 
layers.  From  Figures  4.1  to  4.3,  one  can  realize  that 
the  amount  of  deflection  suddenly  increases  at  the 
location  of  slip  zone.  This  phenomenon  is  pronounced  at 
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the  tiltmeter  furthest  downslope  (BH7).  However,  the 
slope  indicator  data  are  available  for  only  one  year. 

2.  the  sub-surface  ground  movements  are  uncertain. 

The  information  obtained  from  the  field  has  one  major 
disadvantage  and  that  is  that  the  movements  of  these 
stakes  only  represents  the  movement  of  the  ground 
surface . 

3.  regional  ground  water  flow  has  not  been  studied  in 
deta i 1 . 

Local  hydrology  has  not  been  studied  in  detail  for  the 
Edgerton  area.  The  relationship  between  the  change  of 
moisture  content  from  rainfall  and  snowfall  and  the 
change  of  pore  pressure  has  not  been  established  over  a 
long  time  interval. 

In  addition  to  problems  associated  with  the  limited 
field  data,  there  is  a  problem  in  the  numerical  formulation. 
The  jointing  of  the  stiff-fissured  clay  is  similar  to  the 
discontinuity  zones  of  rock.  However,  modelling  the  joint 
set  is  rather  a  difficult  task. 

4.5  Analysis 

The  entire  set-up  for  the  analysis  is  based  on  the  data 
from  field  work  coupled  with  the  suggested  mode  of  failure. 
The  load-transfer  technique  discussed  in  Chapter  3  will  be 
used  to  analyze  the  Edgerton  Slide. 


68 


4.5.1  Mesh 

The  finite  element  mesh  of  the  profile  is  shown  in 
Figure  4.4.  The  positions  of  the  nodes  are  governed  by  the 
locations  of  the  slope  indicators,  the  survey  hubs  and  the 
stratigraphic  profile  of  the  slide.  The  use  of  a  coarse  mesh 
and  of  4-node  elements  was  dictated  by  the  cost  of  computer 
runs . 

Two  thin  layers  of  elements  are  used  to  model  the  shear 
zone  material  and  to  avoid  the  problem  of  a  stress-free 
boundary  (refer  to  Section  3.2.2).  The  element  aspect  ratio 
(which  is  defined  as  the  length  divided  by  the  height  of  the 
element)  ranges  from  0.15  to  500. 

The  boundary  effect  is  minimized  by  setting  the 
boundaries  according  to  the  boundary  studies  by  Desai  and 
Christian  (1977).  The  upslope  boundary  is  approximately  300 
meters  away  from  the  crest.  The  bottom  boundary  is 
approximately  100  meters  below  the  shear  zone.  The  downslope 
boundary  is  approximately  200  meters  away  from  the  central 
probable  rupture  surface. 

Three  probable  rupture  surfaces  as  shown  in  Figure  4.4 
are  considered  in  the  design  of  the  mesh. 

4.5.2  Material  Properties 

Typical  stratigraphy  of  the  site  is  shown  in  Figure 
2.3.  All  of  the  material  parameters  are  based  on  either 
available  data  or  on  a  reasonable  estimate  based  on  local 


experience . 
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Input  parameters  for  numerical  analyses  depend  on  the 
type  of  analyses.  However,  one  of  the  basic  parameters,  mass 
density,  is  required  for  all  analyses.  It  is  used  in  the 
calculation  of  the  element  mass  matrix.  Additional  input 
parameters,  for  example,  are  Young's  Modulus  and  Poisson's 
Ratio.  These  parameters  are  used  to  define  the  material 
properties  for  an  isotropic  linear  elastic  analysis. 

At  the  outset,  all  the  values  of  Young's  Modulus  and 
Poisson's  Ratio  were  derived  from  Balanko  (1981)  and  are 
shown  in  Table  3.3.  The  major  reason  is  that  both  areas 
comprise  similar  material  with  similar  properties  and  a 
similar  stress  history. 

4.5.3  Approach 

For  the  following  numerical  analyses,  both  the 
stratigraphic  and  piezometric  profile  were  assumed  to  be 
unchanged  from  1975. 

The  ratio  of  the  width  to  the  height  of  the  Edgerton 

a 

Slide  is  large,  hence  both  the  side  and  the  end  effects  will 
be  small.  The  results  from  the  two-dimensional  analysis  will 
be  adequate.  Therefore,  the  two-dimensional  analysis  can  be 
used  to  save  computing  time  as  well  as  provide  a  realistic, 
practical  solution. 

The  time-dependent  movement  of  the  Edgerton  landslide 
will  be  analysed  by  using  a  pseudo-elastic  finite  element 
model.  It  is  assumed  in  the  analysis  that  the  Edgerton 
landslide  movement  is  due  to  the  shear  strength  reduction  at 
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the  shear  zone  with  time.  By  varying  the  shear  modulus  (G) 
in  terms  of  the  Young's  Modulus,  different  surface 
displacements  or  horizontal  displacement  at  the  locations  of 
slope  indicator  can  be  obtained  from  the  finite  element 
results.  The  predicted  displacements  are  compared  with  the 
observed  displacements  in  the  field.  The  analytical  approach 
has  been  applied  to  a  simple  model,  as  discussed  in  Chapter 
3. 

This  analysis  assumes  that  the  surface  displacement  is 
predominantly  caused  by  slippage  at  the  shear  zone.  The 
properties  of  the  other  material  (i.e.  except  shear  zone 
material  property)  are  assumed  to  be  t ime- independent . 

4.5.4  Other  Details  For  The  Analysis 

The  original  material  properties  along  the  slip  surface 
elements  are  identical  to  those  of  the  adjacent  elements. 
Afterwards,  the  material  properties  along  the  slip  surface 
are  reduced  uniformly  to  a  lower  value.  The  central  probable 
rupture  surface  as  shown  on  Figure  4.4  is  used  for  joining 
the  interbedded  shear  zone  to  the  surface.  This  surface  is 
favored  over  the  other  two  by  the  field  evidence  of  the 
cracking  pattern  and  the  uplift  movements. 

The  results  from  an  approximate  50  percent  of  Young's 
Modulus  reduction  do  not  agree  in  the  order  of  magnitude  of 
the  field  data  (primarily  displacements).  Therefore,  large 
reduction  of  Young's  Modulus  will  be  used.  The  following 
results  will  be  derived  from  the  reduction  range  of  10  to 
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0.5  percent  of  its  original  value. 

4.6  Results  Of  The  Analyses 

The  horizontal  displacements  from  the  analytical 
results  along  the  locations  of  boreholes  2,  4  and  7  are 
plotted  and  are  shown  in  Figures  4.1  to  4.3.  The  following 
points  are  of  note: 

a.  Changes  in  Young's  Modulus  of  slip  surface  have  a 
reverse  effect  on  the  horizontal  displacements. 

b.  From  a  comparison  of  Figures  4.1,  4.2  and  4.3,  the 
largest  horizontal  displacement  takes  place  in 
borehole  number  2,  which  is  the  furthest  upslope 

t i ltmeter . 

c.  The  areas  within  and/or  adjacent  to  the  shear  zone 
(or  slip  surface)  show  larger  movement  than  those 
below  or  above  it.  This  phenomenon  is  particularly 
evident  in  borehole  number  7,  (Figure  4.3)  which  is 
the  furthest  downslope  tiltmeter. 

Three  maximum  shear  stress  contours  are  shown  in 
Figures  4.5,  4.6  and  4.7.  These  contours  represent  three 
different  stages;  namely,  prior  to  valley  development,  after 
valley  development  and  after  the  development  of  the 
weakening  zone. 

The  maximum  shear  stress  contours  prior  to  valley 
development  are  predominantly  governed  by  the  soil 
properties  and  the  design  of  the  mesh.  (Figure  4.5)  The 
maximum  shear  stress  contours  after  the  development  of  the 
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weakening  zone  (Figure  4.7)  show  a  stress  concentration 
along  the  slip-surface;  in  particular,  the  portion  furthest 
upslope . 

Figures  4.8,  4.9  and  4.10  show  the  vertical,  horizontal 
and  shear  stress  contours  respectively  at  the  stage  after 
the  development  of  the  weakening  zone. 

The  surface  displacements  are  plotted  in  Figure  4.11. 
The  analytical  result  shows  that  the  surface  movements  are 
primarily  moving  downhill. 

4.7  Evaluation  Of  Young's  Modulus  From  The  Previous 

Laboratory  Results 

Tweedie  (1976)  conducted  two  direct  shear  tests  for  the 
remoulded  bentonitic  clayshale.  The  aim  of  this  section  is 
to  derive  the  Young's  Modulus  from  his  laboratory  results. 

The  procedure  which  was  developed  by  Noonan  and  Nixon 
(1972)  will  be  used  to  determine  the  Young's  Modulus  from 
the  direct  shear  test.  The  following  information  together 
with  the  direct  shear  test  results  are  required  to  determine 
the  Young's  Modulus: 

a.  The  sample  size  is  5.08  centimeter  (2  inches)  square 
and  2.54  centimeter  (1  inch)  thich. 

b.  The  gap  between  the  upper  and  lower  halves  of  the 
shear  box  is  approximately  1  millimeter  (by  turning 
the  screws  one  half  to  three  quarter  of  a  turn). 

If  the  Poisson's  Ratio  is  set  to  0.42,  the  Young's 
Modulus  will  range  from  4.0  MPa  to  5  MPa.  Therefore,  the 


Young's  Modulus  obtained  from  the  laboratory  procedure  can 
be  compared  with  that  obtained  from  the  analytical 
procedure . 

4.8  Comparison  And  Discussion  Of  Results 

Some  of  the  analytical  results  can  be  compared  with  the 
available  field  data.  Additionally,  the  input  parameters  for 
the  analysis  can  be  checked  with  the  laboratory  values. 

These  comparisons  are: 

a.  The  horizontal  displacements  from  the  analytical 
results  along  the  locations  of  boreholes  2,  4  and  7 
can  be  compared  with  slope  indicator  readings. 

b.  The  surface  displacement  from  the  analytical  results 
can  be  compared  with  the  field  measurement. 

c.  The  Young’s  Modulus  which  is  used  for  the  analysis 
to  compare  with  the  field  data  can  be  compared  with 
the  one  obtained  from  the  laboratory  results. 

A  discussion  of  the  stress  contours  is  presented  to 
explain  the  failure  mechanism  of  the  Edgerton  Slide. 

4.8.1  Comparison  With  The  Slope  Indicator  Readings 

The  analytical  results  of  both  boreholes  2  and  7  follow 
a  similar  trend  as  the  slope  indicator  readings.  However, 
the  difference  is  that  different  material  properties  were 
used  to  match  the  field  data  in  different  locations.  The 
Young's  Modulus  of  10  MPa  and  1  MPa  were  used  to  match  with 
the  data  of  boreholes  2  and  7  respectively.  These  are  shown 
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in  Figures  4.4  and  4.6.  Perhaps,  this  is  an  indication  that 
the  strength  properties  along  the  shear  zone  are  not 
un i f orm . 

The  analytical  result  of  borehole  4  is  similar  to  the 
slope  indicator  reading  at  the  shear  zone,  but  no  comparison 
can  be  drawn  from  the  movements  above  the  shear  zone.  This 
may  be  due  to  another  rupture  suface  co-existing  with  the 
central  rupture  suface. 

4.8.2  Comparison  With  The  Surface  Displacement 

It  seems  that  the  proposed  analytical  procedure  does 
not  adequately  model  the  heaving  portion  located  about  the 
centre  of  the  slide  mass.  This  is  indicated  by  a  comparison 
of  the  two  profiles  on  Figure  4.11.  However,  if  the  shear 
modulus  above  the  shear  zone  is  reduced  and  the  bulk  modulus 
is  kept  constant,  then  the  analytical  results  may  be 
comparable  to  the  surface  displacements  from  the  field. 

4.8.3  Comparison  Of  The  Value  Of  Young's  Modulus 

The  Young's  Modulus  (  E  j  )  which  was  used  in  the 
analysis  ranges  from  1  MPa  to  10  MPa  across  the  shear  zone 
(or  slip  surface)  elements.  The  results  from  the  above  range 
of  Ej  approximately  match  the  field  data  (primarily  slope 
indicator  readings). 

On  the  other  hand,  the  Young's  Modulus  (Ej  )  from  the 
laboratory  procedure  (direct  shear  test)  was  derived  from 
the  remoulded  bentonitic  clayshale  samples,  which  were  taken 
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from  the  core  of  borehole  number  4  at  the  location  of  the 
failure  plane.  The  Young’s  Modulus  derived  from  the 
laboratory  procedure  is  about  4  to  5  MPa. 

Therefore,  the  Young's  Modulus  of  the  slip  surface 
material  can  be  derived  from  either  the  analytical  approach 
or  the  laboratory  approach.  It  is  because  Ej  falls  in  the 
range  of  E, . 

4.8.4  Discussion  Of  The  Stress  Contours 

It  is  important  to  understand  that  the  stress  contours 
developed  after  the  valley  development  are  represented  by  a 
number  of  horizontal  lines.  This  is  due  to  the  material 
properties  being  uniform  within  each  layer. 

A  rapid  stress  change  can  be  observed  from  the  maximum 
shear  stress  contours  after  the  development  of  the  weakening 
zone  (Figure  4.7);  especially  along  the  up-slope  portion  of 
the  slip  surface. 

The  interpretation  from  the  vertical  stress  contours 
after  the  development  of  the  weaking  zone  (Figure  4.8)  is 
that  a  vertical  stress  transfer  may  be  occurring  along  the 
up-slope  portion  of  the  slip  surface. 

It  seems  that  there  is  little  effect  on  the  horizontal 
stresses  after  the  development  of  the  weakening  zone.  This 
indication  arises  from  the  fact  that  the  horizontal  stress 
contours  (Figure  4.9)  are  essentially  a  number  of  horizontal 
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The  shear  stress  (rxy)  contours  after  the  development 
of  the  weakening  zone  (Figure  4.10)  indicate  that  a  rotation 
of  principal  stresses  may  take  place  along  the  toe  rupture 
surface  and  the  up-slope  portion  of  the  slip  surface.  Note 
that  the  irregularity  of  the  stress  contours  is  primarily 
due  to  the  discretization  of  the  element  stresses. 

4.9  Remarks  And  Summary 

If  the  slope  indicator  readings  were  available  for  more 
than  one  year,  a  further  analytical  reduction  of  the  shear 
strength  can  be  carried  out  in  order  to  compare  the 
analytical  results  with  the  field  measurements  observed  over 
the  longer  time  period. 

The  assumption  of  the  shear  modulus  reduction  at  the 
shear  zone  with  time  seems  to  be  appropriate  in  explaining 
the  failure  mechanism  of  the  Edgerton  Slide.  From  the 
results  of  Figures  4.7  to  4.10,  it  seems  that  the  stress 
concentration  is  higher  in  the  up-slope  part  of  the  slope 
than  the  down-slope  portion.  Therefore,  a  relatively  large 
movement  may  first  take  place  in  the  up-slope  area  of  the 
slide  mass.  Hence,  it  is  possible  that  movement  is 
progressing  from  crest  to  toe. 

One  interesting  conclusion  may  be  drawn  from  this 
analysis.  Since  the  Young's  Modulus  obtained  from  the 
laboratory  procedure  agrees  with  that  obtained  from  the 
analytical  approach,  the  laboratory  value  can  be  used  as  an 
input  parameter  in  the  analytical  analysis  to  reduce  the 
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number  of  trials  to  predict  the  observed  movement. 

Therefore,  the  laboratory  value  can  be  used  in  the 
analytical  analysis  during  the  design  phase  to  predict 
future  movement. 

4.10  Area  Required  For  Further  Research 

Improvements  for  representing  and  analyzing  models  more 
effectively  than  are  already  being  analysed,  are  relatively 
important  and  are  in  demand  by  many  practicing  engineers. 
However,  from  a  research  point  of  view,  the  development  of 
techniques  for  modelling  new  phenomena  , such  as  anisotropy, 
strain  softening  and  dilatancy,  is  rather  important  and 
necessary.  The  formulation  of  the  constitutive  relations  for 
most  geological  materials  encounters  difficulties  and 
requires  a  great  deal  of  rationalization.  Yet,  it  appears  to 
this  writer  that  most  researchers  cannot  distinguish  which 
material  model  is  the  best  for  geological  materials. 

During  the  past  decade,  the  development  of  non-linear 
finite  element  analytical  techniques  seemed  to  be  very 
active.  A  result  is  the  development  of  the  computer  programs 
such  as  ADINA  and  ADINAT.  However,  the  development  of 
non-linear  finite  element  techniques  requires  research  in 
various  areas,  as  mentioned  by  Bathe  (1980),  approximation 
theory,  numerical  methods  and  computer  program 
implementation.  Because  of  this,  the  product  of  this  general 
research  derives  from  researchers  specializing  in  different 
fields.  Therefore,  it  is  important  to  prove  that  the  program 
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does  fullfil  the  original  objectives. 

It  is  the  writer's  opinion  that  the  re 
verification  and  qualification  of  any  non-1 
program  is  as  important  as  the  formulation 
After  the  clarification  of  these  programs, 
research  may  progress  to  another  level. 
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Figure  4.1  Horizontal  Displacement  Of  Borehole  2 
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Figure  4.2  Horizontal  Displacement  Of  Borehole  4 
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Figure  4.3  Horizontal  Displacement  Of  Borehole  7 


FIGURE  4.4  PROFILE  RND  MESH  OF  THE  EDGERTON  74-SOUTH 
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Figure  4.5  Maximum  Shear  Stress  Contours  Of  The  Edgerton 
Slide  (Prior  To  Valley  Development) 
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Figure  4.6  Maximum  Shear  Stress  Contours  Of  The  Edgerton 
Slide  (After  The  Valley  Development) 
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5.  GENERAL  APPLICATIONS  OF  THE  LOAD-TRANSFER  TECHNIQUE  AND 

CONCLUSIONS 

5.1  General  Applications  Of  The  Load-Transfer  Technique 

The  load-t ransf er  technique  was  developed  to  handle  the 
strain-softening  behaviour  associated  with  soils  which  are 
vulnerable  to  progressive  failure.  The  load-transfer  program 
was  shown  in  the  working  stage  and  was  presented  in  Chapter 
3.  This  approach  was  applied  to  the  study  of  the  Edgerton 
Slide  and  was  presented  in  Chapter  4. 

This  program  however  can  be  used  to  solve  other 
problems  associated  with  strain-stiffening  materials  and 
materials  with  time-dependent  effects.  The  concept  can  be 
used  for  stress  analysis  in  a  no-tension  material  as  was 
first  used  by  Zienkiewicz  et.al.  (1968). 

The  load-transfer  technique  is  a  pesudo-elast ic 
analysis.  The  basic  approach  is  to  handle  the  excess  shear 
stress  due  to  whatever  reason  and  the  change  in  elastic 
properties  (  eg.  Young's  Modulus  ( E ) ,  Poisson's  Ratio  (m), 
Shear  Modulus  (G) ,  and  Bulk  Modulus  (K)  )  due  to  whatever 
reason . 

5.1.1  The  Strain-Stiffening  Approach 

This  approach  is  direct  opposite  to  the 
strain-weakening  approach.  However,  the  difference  is  that 
the  former  approach  is  used  for  increasing  shear  modulus  of 
a  material  which  can  sustain  excessive  shear  stresses.  A 
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schematic  diagram  is  used  to  show  the  difference  between 
these  approaches  and  is  shown  in  Figure  5.1.  The 
strain-hardening  approach  may  be  applied  to  a  study  of  the 
swelling  behaviour  of  soils. 

5.1.2  The  Study  Of  Time-dependent  effects 

This  approach  is  a  pseudo-time  dependent  analysis.  This 
can  be  illustrated  by  the  following  example. 

If  the  slope  indicator  readings  of  the  Edgerton  Slide 
were  available  for  more  than  one  year,  then  a  creep  analysis 
could  be  done.  The  shear  modulus  (G)  ,  used  in  calculating 
the  displacements  which  is  comparable  to  the  field 
measurements  of  year  1 ,  can  be  used  to  calculate  the  strain 
at  year  1.  Another  reduction  of  shear  modulus,  due  to 
whatever  reasons,  is  done  in  order  to  compare  the  analytical 
results  with  the  field  measurements  of  year  2.  Finally,  the 
strain  at  year  2  is  calculated.  The  process  of  calculation 
could  be  carried  on  for  year  3  and  so  on.  Ultimately,  a 
strain-time  curve  (creep  curve)  can  be  plotted.  Therefore, 
the  calculated  creep  curve  can  be  classified  as  one  of  the 
following  stages  : 

a.  primary  creep 

b.  steady-state  creep 

c.  tertiary  creep 
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5.1.3  The  Study  Of  No-Tension  Materials 

This  approach  can  be  used  to  study  a  rock  mass  in  its 
natural  state  because  it  usually  cannot  sustain  tension  due 
to  the  presence  of  cracks  and  fissures. 

However,  the  load-transfer  program  has  to  be  modified 
such  that  the  final  stress  (oj)  (refer  to  Appendix  A)  is 
artificially  reduced  to  zero.  This  technique  actually  forces 
the  major  principle  stress  (sign  convention  is  that 
compression  is  negative)  to  zero,  however  the  minor 
principal  stress  and  the  principal  plane  direction  remain 
unchanged . 

5.2  Summary  Of  This  Thesis 

This  research  is  a  documented  case  history  in  applying 
the  analytical  method  to  a  field  problem.  The  prime 
objective  is  to  develop  a  procedure  to  analyze  a  natural 
slope . 

The  advantage  of  numerical  analysis,  as  mentioned  by 
Gibson  (1974),  is  that  this  analysis  can  help  to  distinguish 
among  those  factors  that  are  of  primary  significance  and 
those  that  are  of  secondary  importance. 

Conclusions  from  this  research  concern  four  major 
issues : 

a.  The  correct  usage  of  any  existing  computer  program 
is  emphasized  so  that  meaningful  results  can  be 
obtained.  Sometimes,  the  users  have  to  run  a  simple 
test  on  the  available  function  o‘f  any  program  in 
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order  to  determine  whether  the  program  achieves  this 
goal  or  not. 

b.  The  load-transfer  approach  is  developed  to  model  the 
strain-weakening  material,  which  is  vulnerable  to 
progressive  failure. 

c.  The  achievement  arising  from  this  research  is  that 
the  variation  of  the  Young’s  Modulus  was  used  for 
matching  the  displacement  history  of  the  Slide  mass. 
The  analytical  results  show  that  the  final  results 
were  independent  of  Poisson's  Ratio  and  the  initial 
Young’s  Modulus. 

d.  The  load-transfer  technique  can  be  applied  to  other 
problems  associated  with  strain-hardening  material, 
creep  material  and  no-tension  material. 

5.3  Areas  For  Future  Research 

A  precise  correspondence  between  predicted  and  the 
observed  field  measurements  is  rare.  However,  theoretical 
solutions  aid  in  visualising  possible  failure  mechanisms  in 
different  situations  and  in  developing  sound  judgement 
concerning  stability  problems.  The  following  areas  require 
future  research  : 

a.  A  more  realistic  stress-strain  relationship  should 
be  developed  in  considering  problems  associated  with 
the  dilatancy,  rupture  and  cracking,  and  the  true 
strain-softening  behaviour. 

b.  The  load-transfer  technique  could  be  used  to  study 
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other  similar  slides.  From  the  results  of  the 
numerous  case  studies,  a  general  approach  can  be 
established  for  the  study  of  a  natural  slope, 
c.  The  ADINA  program  could  be  modified  so  that  the 

strain-weakening  stress-strain  behaviour  will  be  one 
of  the  available  material  models.  This  can  increase 
the  efficiency  of  solving  a  problem  and  reduce  the 
amount  of  work  in  preparing  the  input  data  file. 
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STRAIN-STIFFENING  APPROACH 


STRAIN-WEAKENING  APPROACH 


Figure  5.1  Schematic  Diagram  Of  The  Strain-Stiffening 
Approach  And  The  Strain-Weakening  Approach 
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Appendix  A 

The  Finite  Element  Formulation  Of  The 
Incremental  Loading  Due  To 
The  Reduction  Of  Shear  Modulus 
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Appendix  A 

Finite  Element  Formulation  Of  The  Incremental  Loading  Due 
To  The  Reduction  Of  Shear  Modulus 


AT  EQUILIBRIUM 


At  time  equals  T  j , 


the  incremental  finite  element  equilibrium  equation  can  be 
expressed  in  the  following  form  by  using  Virtual 
Displacement  Principle  : 

f[B|]  {o,}dv  =  {Ri} 


where 


N- 
H  - 

{R,}- 


strain  displacement  matrix  at 
internal  stresses  at  time  Tj 
external  load  at  time  Tj 


time  Ti 
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At  time  equals  T i ,  where  Tj  =  Tj  +  6T  , 


due  to  the  change  in  elastic  parameters  at  time  Tj ,  there 
will  be  a  corresponding  changes  in  stress. 

The  strain  (e)  at  time  Tj  is  the  same  as  that  at  the 
beginning  of  time  Tj. 

Therefore,  the  amount  of  stress  change  is  given  by  : 


The  equivalent  nodal  force  due  to  this  change  in 
stresses  is  given  by  : 


This  represent  the  portion  of  the  stress  that  cannot  be 
taken  by  the  weaken  zone.  Hence,  this  stress  must  be 
redistributed  to  the  other  parts  of  the  structure. 
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The  new  equilibrium  equation  becomes  : 


Therefore , 


this  equilibrium  equation  must  be  satisfied  at  the  new 
stress  state. 


Appendix  B 

User's  Manual  and  Source  Code 

Of 

The  Load-Transfer  Program 
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Appendix  B 

USER’S  MANUAL  OF  LOAD-TRANSFER  PROGRAM 

I  n t roduct i on 

The  load-t ransf er  program  consists  of  two  parts;  these 
are  the  main  program  and  the  library  program.  The  following 
MTS  commands  can  be  used  to  run  the  load-transfer  program  : 


$RUN  CLOAD+CLIB  2=IN2  4=IN4  5=IN5  6=-OUT6  7=-OUT7 


where 

CLOAD  =  compiled  version  of  the  main  program 
CLIB  =  compiled  version  of  the  library  program 
IN2  =  input  file  containing  the  modifying  material 
parameters 

IN4  =  file  containing  the  ADINA  output;  primarily 
displacements 

IN5  =  file  containing  ADINA  input 

-OUT6  =  temporary  output  file  echoing  input  information 
-OUT7  =  temporary  output  file  of  the  new  element 
stresses  and  the  redistribution  loads 
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Detail  Descriptions 

Data  For  IN2 

Data  consists  of  : 

1.  FORMAT(I5)  —  the  number  of  elements  to  have  changes 
made  to  their  material  properties. 

2.  FORMAT ( I  5 , 2G 1 0 . 0 )  —  a  list  of  each  element  number, 
new  Young's  Modulus  and  new  Poisson's  Ratio. 

Data  For  IN4 

Data  consi sts  of : 

1.  FORMAT (I7,31X,2G18.6)  -  the  last  step  of  the 
displacement  output  from  ADINA. 

Data  For  IN5 

The  data  format  as  described  in  the  ADINA  manual. 
However,  the  subroutine  INDATA  may  have  to  be  changed 
for  each  material  model  used  other  than  linear-elastic. 

Output  Of  -OUT6 

This  file  is  primarily  used  for  sel f -check ing  input 
i n  f  orma t i on . 

Output  Of  -OUT7 

Output  consists  of  two  parts  : 

1.  the  new  element  stresses  are  calculated  and  printed. 

2.  the  redistribution  loads  are  calculated  and  printed. 
These  will  be  re-used  as  the  input  information  for 


ADINA. 


. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 1 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  MAIN  PROGRAM  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IMPLICIT  REAL*4(A-H,0-Z) 

DIMENSION  B0(6, 16) ,FL( 16) , X ( 8 ) , Y ( 8 ) ,X8(3) ,W8(3) . AN ( 9 ) , ANSI  9) , 

1  ANR ( 9 ) , ANT ( 9 ) 

DIMENSION  ICO ( 10, 1000) ,XX( 1000) ,YY( 1000) .STRESS (4, 4, 1000) , 

1  STRAIN (4, 4, 1000) ,UU( 1000) , VV ( 1 000 ) , ELMP ( 2 , 1000) , 

2  EMPNEW ( 2 , 1000) ,P( 1000) , IX (2000) , NEWEL ( 1000) 

DATA  NICO, INI , IN2, IN3,NGP,NS,NELMP/ 10,5,4,2,4,4,2/ 

DATA  I  OUT  1 , IOUT2/6 , 7/ 

DATA  W8/0 . 55555555555556D0 , 0 . 88888888888889D0 , 

1  0  55555555555556D0/ 

DATA  X8/-0. 77459666924 148D0.0. DO. 0. 77459666924 148D0/ 

NI  P  =  2 

I F ( N I P . GT . 2 )  GO  TO  30 

X8( 1 )=-0. 577350269 1896257D0 

X  8  (  2  )  =  -  X  8  (  1  ) 

W8 ( 1 ) = 1 . DO 
W8 ( 2 ) = 1 . DO 

30  CALL  INDATAdCO.XX, YY.NEL.NNOD, NICO, INI , ELMP, NELMP, 

1  I  OUT  1 , IX) 

CALL  I NSTR ( ICO , XX , YY , NEL , NNOD , STRESS , STRAIN , UU , VV , NS , NGP , 

1  NICO, IN2, IOUT1 ) 

CALL  MODI F( NEWEL, NNEWEL, EMPNEW, IN3, IOUT1 , NELMP) 

CALL  DLOAD ( ICO , BO , FL , X , Y , X8 , W8 , AN , ANS , ANR , ANT , IOUT 1 , NEL , 

1  NNOD, NI P, UU, VV, ELMP, P, EMPNEW, NEWEL, NNEWEL, THICK, NICO, NELMP, 

2  XX  YY  ) 

CALL  RESULT (NNOD.P, IOUT2) 

STOP 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  READ  ADINA  MASTER  INPUT  FILE  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  INDATA ( I  CO , XX , Y Y , NEL , NNOD , N I  CO , I N 1 , ELMP, NELMP, 

1  I  OUT  1  ,  I  X 


YY ( 1 ) , ELMP ( NELMP , 1  )  , I  X ( 1 


IMPLICIT  REAL*4(A-H,0-Z) 

DIMENSION  I  CO ( N I  CO , 1 ) , XX ( 

READ ( INI , 101  )  A 

RE AD ( INI , 102)  NNOD , NGRP 

CALL  I  SET ( IX,2*NNOD) 

WRITE! IOUT1 ,201 )  NNOD, NGRP 
N I  =  1  3 
NN0DEL=8 
DO  1  1  =  1  ,  N I 

1  READ ( INI , 101  )  A 
WRITE ( IOUT 1,202) 

DO  2  1=1 , NNOD 

READ! INI , 103)  I  1 , 12 , XX ( I ) , YY (  I  ) 

I F ( I  1 . EQ . 0 )  I  X ( 2* I  -  1  )  =  1 
I F ( 1 2 . EQ . 0 )  I  X ( 2*  1  )  =  1 

2  WRITE!  IOUT  1  ,  203  )  I  ,  XX  (  I  )  ,  Y  Y  (  I  )  ,  I X  (  2*  I  - 1  )  .  I  X  (  2d 
READ! INI , 101 )  A 

NEL  =  0 


I E  L  =  0 

WRITE! IOUT1 ,204) 


121 

122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 
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137 

138 
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141 
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143 

144 

145 

146 

147 

148 

149 
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157 

158 
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161 
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163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 


1  1  1 


U-DISPL.' , 5X , ' V-DISPL' ,// 


C  WRITE( IOUT1 ,203)  I 
C  DO  3  I G P  = 1 , NGP 

C  READ( IN2, 103)  ( STRESS ( K , IGP , I ) , K= 1 , NS ) 

C  3  WRITE ( IOUT 1 ,204 )  ( STRESS ( K , IGP , I ) , K= 1 , NS : 

C  2  CONTINUE 
RETURN 

101  FORMAT ( I7,31X,2G18.6) 

102  FORMAT ( A4 ) 

103  FORMAT ( 18X.4G15.4) 

201  FORMAT ( // , 5X , ' NODE'  ,5X 

202  F0RMAT(5X, 15.5X.2E15.5) 

203  FORMAT(5X,' STRESSES  OF  ELEMENT  ',15) 

204  FORMAT ( 5X , 4E 1 5 . 5 ) 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  READ  ELEMENTS  ASSOCIATED  WITH  CHANGING  YOUNG'S  MODULUS  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  MODI F ( NEWEL , NNEWEL , EMPNEW , IN3 , I  OUT  1 , NELMP ) 

IMPLICIT  REAL*4( A-H.O-Z ) 

DIMENSION  NEWEL ( 1 ) , EMPNEW ( NE LMP , 1 ) 

READ ( IN3 , 1 0  1  )  NNEWEL 
WRI TE ( IOUT 1,201)  NNEWEL 
DO  1  1=1, NNEWEL 
RE AD ( I N3 ,  102)  NEWEL 


101 

102 

201 

202 


1  WRITE! IOUT 1 , 202 ) 
RETURN 
FORMAT (15) 

FORMAT ( 1 5 , 2G 1 0 . 0 
FORMAT ( // , 5X , ' NO 


I ) , (EMPNEW! J, I ) , J  = 1 . NELMP ) 
NEWEL ( I ) , (EMPNEW! J, I ) , J=1 , NELMP 


OF  ELEMENT  WITH  NEW  STIFFNESS  =  '  ,15) 

FORMAT ( / , 5X , ' ELEMENT  ='  , 1 5 , 5X , ' NEW  ELASTIC  MOD.  = '  , E  1 5 . 5 , 

1  5X , '  NEW  POISSON  RATIO  ='  ,E15.5) 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CALCULATION  OF  DELTA  SIGMA  C 

C  C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


SUBROUTINE  DLOAD ( I  CO , BO , F L , 

1  NEL,NNOD,NIP,UU,VV,ELMP,P, 

2  THICK, NICO, NELMP, XX, YY) 
IMPLICIT  REAL*4 (A-H.O-Z) 
DIMENSION  ICO ( NICO , 1 ) ,B0(6, 

1  AN ( 1 ) , ANS ( 1 ) , ANT ( 1 ) . ANR ( 1 ) 

2  UU ( 1 ) ,VV( 1 ) , ELMP ( NELMP , 1 ) , 
CALL  PSET(P,2*NNOD) 

THICK= 1 .DO 

DO  1  IN=1, NNEWEL 
IEL=NEWEL( IN) 


X , Y , X8 , W8 , AN , ANS , ANR , ANT , IOUT 1 
EMPNEW, NEWEL, NNEWEL, 


1  )  ,FL(  1  )  ,X<  1  ) 
,U( 16) ,EMP(2) 
P ( 1 ) , NEWEL ( 1 ) 


Y ( 1 ) , X8 ( 1 ) , W8 ( 1 ) 
EMPN ( 2 ) , XX ( 1 ) , YY 
EMPNEW ( NELMP , 1 ) 


1  ) 


CALL 
CALL 
CALL 
DO  2 


PSET 
PSET 
PSET 
11  =  1 


I  ICO= ICO ( I  1 , I  EL ) 
IF( IICO.EQ.O)  GO 
X( 1 1 )=XX( I  ICO) 

Y ( I  1 ) =YY ( I  ICO ) 
U(2*1 1-  1  )  =UU ( I  ICO 
U ( 2* 1 1 )=VV( I  ICO ) 


TO  2 
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2  CONTINUE 

DO  3  11=1,2 
I IC0= ICO ( 1 0 . 1  EL ) 

EMP ( 1 1 ) =ELMP ( 1 1 , 1  ICO ) 

3  EMPN ( 1 1 )=EMPNEW( 1 1 , IN) 

C  WRITE (6.1201)  ( X ( 1 1  ) , 11  =  1,8) 

1201  FORMAT ( / , 5X , ' X  =  ' .8F10.3) 

C  WRITE(6, 1202)  ( Y ( 1 1 ) , 1 1  =  1 , 8 ) 

1202  FORMAT ( / , 5X , ' Y  =  ' .8F10.3) 

C  WRITE ( 6 , 1203)  ( U ( 1 1 ) , 1 1  =  1 , 1 6 ) 

1203  FORMAT ( / , 5X , ' U  =  ' , 2 ( 8E 12 . 5 , / , 9X ) ) 

C  WRITE ( 6 , 1204 )  ( EMP ( 1 1  )  ,  1 1  =  1  , 2 ) 

1204  FORMAT ( / , 5X , ' EMP  =  '.8F10.3) 

C  WRITE ( 6 , 1205 )  ( EMPN ( 1 1 > , 1 1  =  1 , 2 ) 

1205  FORMAT ( / , 5X , ' EMPN  =  ' , 8F 10 . 3 ) 

CALL  LOAD ( IEL , BO, FL , X , Y , ICO.NEL , X8 , W8, NIP , AN , ANS, ANT , 

1  IOUT1 .THICK, EMP, EMPN, NI CO, ANR.U) 

CALL  SETUP( P, FL , ICO.NICO, IEL ) 

1  CONTINUE 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  WRITING  THE  REDISTRIBUTION  LOADS 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  RESULT ( NNOD , P , IOUT2 ) 

IMPLICIT  REAL*4 ( A- H , O-Z ) 

DIMENSION  P ( 1 ) 

DO  1  I NOD= 1 , NNOD 
11=2*1  NOD- 1 
I2=2*INOD 
K 1  =  1 
K2  =  2 
K3  =  3 
A  1 =0 . 00 

WRITE ( IOUT2 , 201 ) 

1  WRITE! I OUT2 ,201 ) 

RETURN 

201  FORMAT (3I5.2E10.3) 

END 


I  NOD , K2 , K 1 , P ( I  1 ) , A  1 
INOD.K3.K1 ,P( 12) , A  1 


file 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  LIBRARY  FOR  LOAD-TRANSFER  PROGRAM  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CALCULATION  OF  EXTERNAL  LOAD  FOR  EACH  ELEMENT  C 

C  (  R  MATRIX  )  C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  LOAD ( I  EL , BL , FL , X , Y , ICO , NEL , 

1  X8 , W8 , NI P , AN , ANS , ANT , I  OUT  1 , THICK , EMP , EMPN , NICO , ANR , U ) 
IMPLICIT  REAL*4 ( A - H , 0 - Z ) 

DIMENSION  BL ( 6  t 1 ) ,FL( 1 ) ,X( 1 ) , Y( 1 ) ,Z< 1 ) , EMP ( 1 ) ,EMPN( 1 ) , 

1  I CO (NICO, 1 ) ,X8( 1 ) ,W8( 1 ) ,AN( 1 ) ,ANS( 1 ) ,ANT ( 1 ) , 

2  A J (3,3) , A I (3,3) , ANR ( 1 ) , CE (3,3) , STR ( 3 ) ,FO(3) ,F 1 (3) , 

3  SIG ( 6 ) , F ( 16) ,U( 1 ) 

CALL  PRESET(BL,6, 16) 

CALL  PSET ( SIG , 6 ) 

CALL  PSET ( FL , 16) 

I G  P  =  1 


DO  26  1  =  1  , NIP 
R=X8 ( I ) 

DO  26  J= 1 , N I P 
S  =  X8 ( J ) 

CALL  SHAPE (ANR, ANS, R,S, I  EL, I  CO , NEL , 8 , AN , N I  CO , 2 ) 

CALL  BOM AT (AN,ANR,ANS,ANT,X,Y,BL,AJ,AI ,DET, IEL, I  OUT  1 , 

1  8, 2, Z, 2, RAD) 

CALL  MULT3(BL,6,U,STR,3, 16,1) 

C  WRITE* I0UT1 , 1202)  ( STR ( 1 1 ) , 1 1  =  1 , 3 ) , DET 

1202  FORMAT)/, 5X,' STR  =  ' , 3E 1 2 . 5 , 5X , ' DET  =' , E 15 . 5 ) 

CALL  ELASTC(CE , EMP , IEL ) 

CALL  MULT3(CE,3,STR,FO,3,3, 1 ) 

C  WRITE* I0UT1 , 1203)  ( FO ( 1 1 ) , 1 1 = 1 , 3 ) 

1203  FORMAT ( / , 5X , '  FO  =  '  , 3E 12 . 5 ) 

CALL  ELASTCICE, EMPN, IEL) 

CALL  MULT3(CE,3,STR,F1 ,3,3, 1 ) 

C  WRITE ( 6 , 1201 )  ( (BL(I,J) ,1=1,6) ,J=1, 16) , (FO(I ) ,1=1,3) 

1201  FORMAT*/, 5X, ' BL  IN  LOAD  =  ' , 1 6 ( / , 5X , 6E 1 2 . 5 ) , / , 5X , ' SIG  IN  LOAD 

1  , / , 5X , 3E 1 2 . 5 ) 

WRITE* IOUT1 , 1204)  ( F 1 ( 1 1 ) , 1 1 = 1 , 3 ) 

1204  FORMAT ( / , 5X , ' FI  =  ' .3E12.5) 

DO  3  K=1 ,3 

3  SIG ( K ) =FO ( K ) - F 1 ( K ) 

CALL  MULT3(BL,6,SIG,F,3, 16,2) 

DO  4  K=  1  , 1 6 

4  FL(K)=FL(K)+DET*W8( I )*W8(d)*F (K)*THICK 
26  IGP= IGP+ 1 

C  WRITE* IOUT1 ,201 )  ( FL ( I ) , I = 1 , 1 6 ) 

201  FORMAT(/,5X,' EQUIVALENT  LOAD  VECTOR  IS  '  , 2 ( / , 5X , 8E  1 5 . 6 ) ) 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CALCULATION  OF  SHAPE  FUNCTIONS  C 

C  C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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SUBROUTINE  SHAPE ( ANR , ANS , R , S , I  EL , ICO , NEL , NODE , AN , NICO . I  PS ) 
IMPLICIT  REAL*4 ( A-H , 0- Z ) 

DIMENSION  ANR( 1 ) ,ANS( 1 ) , ICO (NICO, 1 ) , AN ( 1 ) 

ANR(9)=-2.D0*R*( 1 .D0-S*S) 

ANS( 9 ) =-2 . D0*S* ( 1 . D0-R*R ) 

AN ( 9 ) = ( 1 . DO- R*R ) * ( 1 .D0-S*S) 

IF (NODE. EQ. 9)  GO  TO  4 
ANR I  9 ) =0 . DO 
ANSI  9 ) =0 . DO 
AN ( 9 ) =0 . DO 

4  ANR I  7 )  =  -  R* ( 1 .D0+S)-ANR(9)/2.D0 
ANR ( 8 ) = - ( 1 .D0-S*S)/2.D0-ANR(9)/2.D0 
ANR(5)=-R*( 1 .D0-S)-ANR(9)/2.D0 
ANR ( 6 ) = ( 1 .D0-S*S)/2.D0-ANR(9)/2.D0 
ANS ( 7 ) = ( 1 .D0-R*R)/2.D0-ANS(9)/2.D0 
ANS(8)=-S*( 1 .DO-R ) -ANSI9  )/2 .DO 
ANS I  5 )= -  I 1 .DO-R*R ) /2. DO-ANSI  9 )/2. DO 
ANS ( 6 ) = - S* ( 1 . DO+R ) -ANS I  9 ) /2 . DO 
AN  I  7 )  =  I  1 ,DO-R*R)*( 1 .D0+S)/2.D0-AN(9)/2.D0 
AN ( 8 )  =  I  1 .DO-S*S)*( 1 .D0-R)/2.D0-AN(9)/2.D0 
AN  I  5 )  =  I  1 . DO- R*R ) *  I 1 .D0-S)/2.D0-AN(9)/2.D0 
AN  I  6 )  =  I  1 .DO-S*S)*( 1 .DO+R  )/2. DO- AN  I  9) /2. DO 
DO  2  1=5,8 
IFI ICO! I , IEL  )  )  1,1,2 

1  ANR ( I ) =0 . DO 
ANSI  I )=O.DO 
AN ( I ) =0 . DO 

2  CONTINUE 

ANR I  3 )  =  I  1 .D0+S)/4.D0-(ANR(6)+ANR(7) ) /2 . DO- ANR I  9 ) /4 . DO 
ANR I  4 ) = -  I  1 .D0+S)/4.D0-(ANR(7)+ANR(8) ) /2 . DO* ANR ( 9 ) /4 . DO 
ANRI  1 )  =  -( 1 .D0-S)/4.D0- (ANR(5)+ANR(8) ) /2 . DO- ANR ( 9 ) /4 . DO 
ANR I  2 )  =  I  1 . DO- S ) /4 . DO- I ANR ( 5  ) +ANR I  6 ) ) /2 . DO- ANR ( 9 ) /4 . DO 
ANS I  3 )  =  I  1 .D0+R)/4.D0-(ANS(6)+ANS(7) ) /2 . DO- ANS I  9 ) /4 . DO 
ANS I  4 )  =  ( 1 .D0-R)/4.D0-(ANS(7)+ANS(8) ) /2 . DO-ANSI  9 ) /4 . DO 
ANSI  1 )  =  -( 1 .D0-RI/4.D0- I ANS I  5 ) +ANS I  8 ) ) /2 . DO- ANS I  9 ) /4 . DO 
ANSI  2 ) = - I  1 .D0+R)/4.D0- I ANS I  5 ) +ANS I  6 ) ) /2 . DO- ANS I  9 ) /4 . DO 
AN  I  3 )  =  I  1 .DO+R)*( 1 .D0+S)/4.D0-(AN(6)+AN(7) ) /2 . DO- AN  I  9 ) /4 . DO 
AN  I  4 )  =  ( 1 . DO- R ) *  I  1 .D0+S)/4.D0-(AN(7)+AN(8)  ) /2 . DO- AN  I  9 ) /4 . DO 
AN! 1 )  =  ( 1 . DO-R ) *  I  1 . DO- S ) /4 . DO- I  AN  I  5 ) +AN I  8 ) ) / 2 . D 0 - AN  I  9 ) / 4  .  D 0 
AN  I  2 )  =  I  1 . DO  +  R ) *  I  1 .D0-S)/4.D0-(AN(5)+AN(6) ) /2 . DO- AN  I  9 ) /4 . DO 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CALCULATION  OF  B  MATRIX  ---  STRAIN-DISPLACEMENT 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  BOMAT I  AN , ANR . ANS , ANT , X , Y , BL , A J , A I , DET , I  EL , I  OUT  1 , 

1  NODE, IPS, Z.NVAR.R) 

IMPLICIT  REAL*4 1  A - H , 0 - Z ) 

DIMENSION  ANSI  1 ) ,ANT(  1  )  ,X(  1  ),Y( 1 ) ,BL(6, 1 ) ,AJ(3, 1 ) ,AI(3, 1 ) . 

1  Z I  1 ) , ANR I  1  )  , AN  I  1  ) 

CALL  PRESET (AJ.3,3) 

CALL  PRESET ( AI ,3,3) 

R=0 . DO 

DO  1  K  = 1 , NODE 

IFI IPS. EQ. 3)  R=R+AN I K ) *X I K ) 

A J I  1 , 1 ) = A J I  1 , 1 ) +ANR I K ) *X I K  ) 

A J I  1 , 2 ) = A J I  1 , 2 ) +ANR I K ) *  Y I K ) 

A J I  2 , 1  )  = A J I  2  ,  1 ) +  ANS  I K ) *X I K  ) 
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AJ(2,2)=AJ(2,2)+ANS(K)*Y(K) 

1  CONTINUE 

3  DE T  =  A J ( 1 , 1 )*Ad(2,2)-Ad( 1 ,2)*Ad(2, 1 ) 
IF ( DET . LE . 0 . DO  )  GO  TO  999 

AI(1, 1 ) = A J (2.2) /DET 
A I ( 1 , 2 )  =  - Ad ( 1,2) /DET 
AI (2, 1 )  =  -  Ad ( 2 . 1 ) /DET 
A I ( 2 , 2  )  = Ad ( 1 , 1  ) /DET 

4  DO  2  K=1 .NODE 

K 1 =NVAR* ( K- 1 )+1 


DX=AI (1,1 )*ANR(K)+AI ( 1 , 2 ) *ANS ( K ) +AI ( 1 ,3)*ANT(K) 

DY=AI (2, 1 )*ANR(K)+AI(2,2)*ANS(K)+AI(2,3)*ANT(K) 

BL ( 1 , K 1 ) =DX 
BL ( 2 , K 1+ 1 ) =DY 
BL ( 3 , K 1 )=DY 
BL ( 3 , K 1 + 1 )=DX 
IF ( I  PS-3 )  2.5,6 

5  BL ( 4 , K  1  ) =  AN ( K  )  / R 
GO  TO  2 

6  DZ=AI (3,1 )*ANR(K)+AI(3,2)*ANS(K)+AI(3,3)*ANT(K) 

BL ( 4 , K 1 +2  )  =DZ 

BL ( 5 , K 1 + 1 ) =  D  Z 
BL ( 5 , K 1 +2 ) =DY 
BL ( 6 , K 1 ) =DZ 
BL ( 6 , K 1 +2 ) =DX 
2  CONTINUE 
RETURN 

999  WR I TE ( IOUT 1 ,201 )  I  EL 

201  F0RMAT(5X,'  *****ERROR*****  DET.  OF  JACOBIAN  ', 

1  'MATRIX  IS  ZERO  ELEMENT  NUMBER  =  ',15) 

STOP 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CHECKING  THE  NODAL  FORCE  EQUILIBRIUM  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  SE TUP ( B , F L , ICO , N I  CO . I  EL ) 

IMPLICIT  REAL*4( A-H.O-Z) 

DIMENSION  B ( 1 ) , FL ( 1 ) , ICO ( NI CO , 1 ) 

DO  1  1=1,8 
1 1 C0= I  CO ( I , I  EL ) 

IF(IICO.EQ.O)  GO  TO  1 

B(2*I ICO- 1 ) =B(2*IIC0- 1 )+FL (2*1-1 ) 

B(2*IIC0)=B(2*IIC0)+FL(2*I  ) 

1  CONTINUE 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  MATRIX  MULTIPLICATION  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  MULT  3 ( X , NXX , Y , Z , NX , NY , I  CODE ) 

X ( NX , N Y  )  ,  BY  A  VECTOR  Y ( 1  ) 


C  MULTIPLY  A 

C  IF  ICODE= 1 

C  IF  IC0DE=2 

IMPLICIT 


2-D  MATRIX 

X*  y 

Y ( TRANSPOSE ) *X 
REALM  (A-H.O-Z) 


OR  X(TRANSPOSE)*Y 


DIMENSION  X ( NXX ,  1 
IF ( I  CODE .NE .  1  )  GO 


Y  (  1  )  ,  Z  (  1  ) 
TO  3 
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9  DO  1  I X= 1 , NX 
XX=0 . DO 
DO  2  I Y  = 1 , NY 

2  XX=XX+X(IX,IY)*Y( I Y  ) 

1  Z(IX)=XX 
RETURN 

3  DO  4  I X= 1 , NY 
XX=0 . DO 

DO  5  I Y  = 1 , NX 
5  XX=XX+X( IY, IX)*Y( IY) 

4  Z ( I X ) =XX 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  BOTH  PRESET  AND  PSET  ARE  USED  TO  SET  ZERO  MATRIX  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  PRESET ( A , M , N ) 

IMPLICIT  REAL*4( A-H.O-Z) 

DIMENSION  A ( M , 1 ) 

DO  1  1=1 , M 
DO  2  J=  1  ,N 

2  A (  I  ,  J  )  =0 . DO 

1  CONTINUE 
RETURN 
END 

SUBROUTINE  PSET(A.N) 

IMPLICIT  REAL*4(A-H,0-Z) 

DIMENSION  A ( 1 ) 

DO  2  J=1 , N 

2  A ( J ) =0 . DO 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CONSTITUTIVE  MATRIX  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  ELASTC ( CE , EMP , I  EL ) 

•H.O-Z) 

EMP ( 1 ) 

IM.NELMP, ( E  L  M  P ( I , IM) , 1  =  1 , NELMP ) 
IM.NELMP.ELMP' . / ,5X, 315, / .5X.6E12.5) 


1201 


IMPLICIT  REAL*4 ( A- 
DIMENSION  CEO,  1  ) 
WRITE (6, 1201  )  I  EL 
FORMAT ( // , 5X , '  I  EL 
CALL  PRESET ( CE , 3 , 3 ) 
E=EMP ( 1 ) 

V=EMP ( 2 ) 

C  =  E / ( ( 1 . DO  +  V ) ^ ( 1 .DO-2 


DO*  V 


C 1  =  (  1  .DO 
C2= V*C 
( 1 .DO-2 
4  1=1,2 
3  d=  1  ,2 
I , J ) =C2 
I  ,  I  )  =  C  1 
CONTINUE 
CE ( 3 , 3 ) =C3 
RETURN 
END 


V)*C 


C3 

DO 

DO 

CE 

CE 


DO*V ) *C/2 . DO 
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